IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19532 for trunk/ppStack


Ignore:
Timestamp:
Sep 11, 2008, 6:12:42 PM (18 years ago)
Author:
Paul Price
Message:

Correcting stack variances for the imperfect convolution. Changed (again) the way the deviations are calculated; they're now straight chi2, which allows us to estimate how the variance factor needs to be changed in order to correct for it. Moving additional values into pmSubtractionKernels so that this can be calculated from a kernel that has been read in (very useful for debugging, also nice to have a record of the quality of the fit).

Location:
trunk/ppStack/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r19346 r19532  
    8989                               const psArray *readouts, // Input readouts
    9090                               const psArray *regions, // Array with array of regions used in each PSF match
    91                                const psArray *kernels // Array with array of kernels used in each PSF match
     91                               const psArray *kernels, // Array with array of kernels used in each PSF match
     92                               const psVector *addVariance // Additional variance for rejection
    9293    );
    9394
    … …  
    128129
    129130/// Convolve image to match specified seeing
    130 bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output
    131                   psArray **regions, // Array of regions used in each PSF matching, returned
    132                   psArray **kernels, // Array of kernels used in each PSF matching, returned
     131bool ppStackMatch(pmReadout *readout,   // Readout to be convolved; replaced with output
     132                  psArray **regions,    // Array of regions used in each PSF matching, returned
     133                  psArray **kernels,    // Array of kernels used in each PSF matching, returned
     134                  float *chi2,          // Chi^2 from the stamps
    133135                  const psArray *sources, // Array of sources
    134136                  const pmPSF *psf,     // Target PSF

  • trunk/ppStack/src/ppStackLoop.c

    r19505 r19532  
    390390    psVector *inputMask = psVectorAlloc(num, PS_TYPE_U8); // Mask for inputs
    391391    psVectorInit(inputMask, 0);
     392    psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
     393    psVectorInit(matchChi2, NAN);
    392394    for (int i = 0; i < num; i++) {
    393395        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    … …  
    402404            psFree(rng);
    403405            psFree(inputMask);
     406            psFree(matchChi2);
    404407            return false;
    405408        }
    … …  
    419422            psFree(rng);
    420423            psFree(inputMask);
     424            psFree(matchChi2);
    421425            return false;
    422426        }
    … …  
    425429        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    426430        psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching
    427         if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, rng, config)) {
     431        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i],
     432                          sources, targetPSF, rng, config)) {
    428433            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    429434            inputMask->data.U8[i] = PPSTACK_MASK_MATCH;
    … …  
    459464        psFree(cells);
    460465        psFree(inputMask);
     466        psFree(inputMask);
     467        psFree(matchChi2);
    461468        return false;
    462469    }
    … …  
    473480    memDump("preinitial");
    474481
    475     // Stack the convolvd files
     482    // Stack the convolved files
    476483    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
    477484    pmReadout *outRO = NULL;            // Output readout
    … …  
    489496            psFree(stack);
    490497            psFree(inputMask);
     498            psFree(matchChi2);
    491499            psFree(exptimes);
    492500            psFree(cells);
    … …  
    502510            psFree(stack);
    503511            psFree(inputMask);
     512            psFree(matchChi2);
    504513            psFree(exptimes);
    505514            psFree(cells);
    … …  
    518527            psFree(stack);
    519528            psFree(inputMask);
     529            psFree(matchChi2);
    520530            psFree(view);
    521531            psFree(outRO);
    … …  
    546556        psFree(cellList);
    547557
    548     psFree(cells);
     558        psFree(cells);
    549559
    550560        bool status;                    // Status of read
    … …  
    559569                psFree(stack);
    560570                psFree(inputMask);
     571                psFree(matchChi2);
    561572                psFree(view);
    562573                psFree(outRO);
    … …  
    576587            psArrayAdd(job->args, 1, subRegions);
    577588            psArrayAdd(job->args, 1, subKernels);
     589            psArrayAdd(job->args, 1, matchChi2);
    578590            if (!psThreadJobAddPending(job)) {
    579591                psFree(job);
    … …  
    582594                psFree(stack);
    583595                psFree(inputMask);
     596                psFree(matchChi2);
    584597                psFree(view);
    585598                psFree(outRO);
    … …  
    596609            psFree(stack);
    597610            psFree(inputMask);
     611            psFree(matchChi2);
    598612            psFree(view);
    599613            psFree(outRO);
    … …  
    601615            return false;
    602616        }
     617        psFree(matchChi2);
    603618
    604619        // Harvest the jobs, gathering the inspection lists

  • trunk/ppStack/src/ppStackMatch.c

    r19486 r19532  
    1515                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
    1616
    17 //#define TESTING                         // Enable debugging output
     17#define TESTING                         // Enable debugging output
    1818
    1919
    … …  
    4848#endif
    4949
    50 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
     50bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2,
    5151                  const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
    5252{
    … …  
    111111                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    112112            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
    113             psMetadataItem *vfItem = psMetadataLookup(output->parent->concepts, "CELL.VARFACTOR");
     113            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
    114114            if (!isfinite(vf)) {
    115115                vf = 1.0;
    … …  
    341341    assert((*regions)->n == (*kernels)->n);
    342342
     343    // Correct the variance for the chi^2
     344    {
     345        *chi2 = 0.0;
     346        int num = 0;                    // Number of measurements of chi^2
     347        psString regex = NULL;          // Regular expression
     348        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     349        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
     350        psFree(regex);
     351        psMetadataItem *item = NULL;// Item from iteration
     352        while ((item = psMetadataGetAndIncrement(iter))) {
     353            assert(item->type == PS_DATA_UNKNOWN);
     354            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
     355            *chi2 += kernels->mean;
     356            num++;
     357        }
     358        psFree(iter);
     359
     360        float vf = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.VARFACTOR"); // Variance factor
     361        *chi2 /= vf * num;
     362    }
     363
    343364    // Renormalise the variances if desired
    344365    if (renorm) {

  • trunk/ppStack/src/ppStackReadout.c

    r19475 r19532  
    2222    psArray *subRegions = args->data[3]; // Regions for PSF-matching
    2323    psArray *subKernels = args->data[4]; // Kernels for PSF-matching
    24 
    25     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels);
     24    psVector *addVariance = args->data[5]; // Additional variance when rejecting
     25
     26    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts,
     27                                             subRegions, subKernels, addVariance);
    2628
    2729    job->results = inspect;
    … …  
    7577
    7678psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    77                                const psArray *regions, const psArray *kernels)
     79                               const psArray *regions, const psArray *kernels, const psVector *addVariance)
    7880{
    7981    assert(config);
    … …  
    8486    assert(readouts->n == regions->n);
    8587    assert(regions->n == kernels->n);
     88    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
    8689    static int sectionNum = 0;          // Section number; for debugging outputs
    8790
    … …  
    127130        }
    128131
    129         stack->data[i] = pmStackDataAlloc(ro, weighting);
     132        stack->data[i] = pmStackDataAlloc(ro, weighting, addVariance->data.F32[i]);
    130133    }
    131134
    … …  
    221224        }
    222225
    223         pmStackData *data = pmStackDataAlloc(ro, weighting);
     226        pmStackData *data = pmStackDataAlloc(ro, weighting, NAN);
    224227        data->reject = psMemIncrRefCounter(rejected->data[i]);
    225228        stack->data[i] = data;

  • trunk/ppStack/src/ppStackThread.c

    r19346 r19532  
    235235
    236236    {
    237         psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 5);
     237        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6);
    238238        task->function = &ppStackReadoutInitialThread;
    239239        psThreadTaskAdd(task);
Note: See TracChangeset for help on using the changeset viewer.