IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20711


Ignore:
Timestamp:
Nov 12, 2008, 4:59:10 PM (17 years ago)
Author:
Paul Price
Message:

Can't figure out why the weightings which are being placed on readout->analysis aren't getting to the combination stage, so I'm putting them in a vector and passing them explicitly.

Location:
trunk/ppStack/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r20659 r20711  
    9999                               const psArray *regions, // Array with array of regions used in each PSF match
    100100                               const psArray *kernels, // Array with array of kernels used in each PSF match
     101                               const psVector *weightings, // Weighting factors for each image
    101102                               const psVector *addVariance // Additional variance for rejection
    102103    );
    … …  
    114115                         pmReadout *outRO,   // Output readout
    115116                         const psArray *readouts, // Input readouts
    116                          const psArray *rejected // Array with pixels rejected in each image
     117                         const psArray *rejected, // Array with pixels rejected in each image
     118                         const psVector *weightings // Weighting factors for each image
    117119    );
    118120
    … …  
    142144                  psArray **kernels,    // Array of kernels used in each PSF matching, returned
    143145                  float *chi2,          // Chi^2 from the stamps
     146                  float *weighting,     // Stack weighting (1/noise^2)
    144147                  const psArray *sources, // Array of sources
    145148                  const pmPSF *psf,     // Target PSF

  • trunk/ppStack/src/ppStackLoop.c

    r20659 r20711  
    444444    psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
    445445    psVectorInit(matchChi2, NAN);
     446    psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
     447    psVectorInit(weightings, NAN);
    446448    for (int i = 0; i < num; i++) {
    447449        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    … …  
    482484        psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching
    483485        psTimerStart("PPSTACK_MATCH");
    484         if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i],
     486        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
    485487                          sources, targetPSF, rng, config)) {
    486488            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    … …  
    718720            psArrayAdd(job->args, 1, subRegions);
    719721            psArrayAdd(job->args, 1, subKernels);
     722            psArrayAdd(job->args, 1, weightings);
    720723            psArrayAdd(job->args, 1, matchChi2);
    721724            if (!psThreadJobAddPending(job)) {
    … …  
    725728                psFree(stack);
    726729                psFree(inputMask);
     730                psFree(weightings);
    727731                psFree(matchChi2);
    728732                psFree(view);
    … …  
    961965            psArrayAdd(job->args, 1, outRO);
    962966            psArrayAdd(job->args, 1, rejected);
     967            psArrayAdd(job->args, 1, weightings);
    963968            if (!psThreadJobAddPending(job)) {
    964969                psFree(job);

  • trunk/ppStack/src/ppStackMatch.c

    r20710 r20711  
    4747#endif
    4848
    49 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2,
     49bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
    5050                  const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
    5151{
    … …  
    384384        (void)psBinaryOp(readout->image, readout->image, "+",
    385385                         psScalarAlloc(- psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
     386        *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
    386387        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
    387388                         "Weighting by 1/noise^2 for stack",

  • trunk/ppStack/src/ppStackReadout.c

    r20710 r20711  
    2222    psArray *subRegions = args->data[3]; // Regions for PSF-matching
    2323    psArray *subKernels = args->data[4]; // Kernels for PSF-matching
    24     psVector *addVariance = args->data[5]; // Additional variance when rejecting
    25 
    26     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts,
    27                                              subRegions, subKernels, addVariance);
     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,
     28                                             weightings, addVariance);
    2829
    2930    job->results = inspect;
    … …  
    4243    pmReadout *outRO = args->data[2];   // Output readout
    4344    psArray *rejected = args->data[3];  // Rejected pixels
    44 
    45     bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected); // Status of operation
     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
    4649
    4750    thread->busy = false;
    … …  
    7780
    7881psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    79                                const psArray *regions, const psArray *kernels, const psVector *addVariance)
     82                               const psArray *regions, const psArray *kernels, const psVector *weightings,
     83                               const psVector *addVariance)
    8084{
    8185    assert(config);
    … …  
    8690    assert(readouts->n == regions->n);
    8791    assert(regions->n == kernels->n);
     92    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    8893    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
    8994    static int sectionNum = 0;          // Section number; for debugging outputs
    … …  
    118123        }
    119124
     125#if 0
     126        // This doesn't seem to work, so getting the weightings directly from a vector
    120127        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
    121128        if (!mdok || !isfinite(weighting)) {
    … …  
    125132            psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting);
    126133        }
     134#endif
    127135
    128136        // Ensure there is a mask, or pmStackCombine will complain
    … …  
    132140        }
    133141
    134         stack->data[i] = pmStackDataAlloc(ro, weighting, addVariance->data.F32[i]);
     142        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
    135143    }
    136144
    … …  
    176184
    177185bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    178                          const psArray *rejected)
     186                         const psArray *rejected, const psVector *weightings)
    179187{
    180188    assert(config);
    … …  
    183191    assert(rejected);
    184192    assert(readouts->n == rejected->n);
     193    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    185194
    186195    static int sectionNum = 0;
    … …  
    212221        assert(ro);
    213222
     223#if 0
     224        // This doesn't seem to work, so getting the weightings directly from a vector
    214225        bool mdok;                      // Status of MD lookup
    215226        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
    … …  
    218229            weighting = 1.0;
    219230        }
     231#endif
    220232
    221233        // Ensure there is a mask, or pmStackCombine will complain
    … …  
    225237        }
    226238
    227         pmStackData *data = pmStackDataAlloc(ro, weighting, NAN);
     239        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN);
    228240        data->reject = psMemIncrRefCounter(rejected->data[i]);
    229241        stack->data[i] = data;

  • trunk/ppStack/src/ppStackThread.c

    r19532 r20711  
    235235
    236236    {
    237         psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6);
     237        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 7);
    238238        task->function = &ppStackReadoutInitialThread;
    239239        psThreadTaskAdd(task);
    … …  
    249249
    250250    {
    251         psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 4);
     251        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 5);
    252252        task->function = &ppStackReadoutFinalThread;
    253253        psThreadTaskAdd(task);
Note: See TracChangeset for help on using the changeset viewer.