IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26007 for branches/pap/ppStack


Ignore:
Timestamp:
Nov 2, 2009, 5:08:31 PM (17 years ago)
Author:
Paul Price
Message:

Serious rework of pmStackCombine for clarity and some additional functionality. Split the combinePixels function into two operating modes (test and combine) to make clear what we're doing; pmStackCombine could be split into two separate functions (e.g., pmStackTest and pmStackCombine) in the future, since the testing part doesn't need to produce any output image, but this involves moving a lot of stuff around in ppStack, so not doing it yet. Verified that variance scaling is correct, and softening of variance is correct. Big improvement in minimising loss of good pixels comes from throwing out only the most variant pixel on each rejection pass (consequently changed iterations parameter to iterations per input); this involved another rework of the core combinePixels function (moving the switch on the number of valid inputs inside the iteration). Much happier with the outlier rejection now.

Location:
branches/pap/ppStack/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/ppStack/src/ppStackArguments.c

    r23841 r26007  
    148148    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
    149149    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
    150     psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
     150    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-iter", 0, "Number of rejection iterations per input", NAN);
    151151    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0,
    152152                     "Combination rejection thresold (sigma)", NAN);
    … …  
    250250    }
    251251
    252     VALUE_ARG_RECIPE_INT("-iter",              "ITER",            S32, 0);
     252    VALUE_ARG_RECIPE_FLOAT("-combine-iter",    "COMBINE.ITER",    F32);
    253253    VALUE_ARG_RECIPE_FLOAT("-combine-rej",     "COMBINE.REJ",     F32);
    254254    VALUE_ARG_RECIPE_FLOAT("-combine-sys",     "COMBINE.SYS",     F32);
    … …  
    260260    VALUE_ARG_RECIPE_FLOAT("-poor-frac",       "POOR.FRACTION",   F32);
    261261
    262     valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);
     262    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.IN",   recipe);
    263263    valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
    264264    valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);

  • branches/pap/ppStack/src/ppStackCombineFinal.c

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING                         // Enable test output
     12//#define TESTING                         // Enable test output
    1313
    1414bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,

  • branches/pap/ppStack/src/ppStackCombineInitial.c

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING                         // Enable test output
     12//#define TESTING                         // Enable test output
    1313
    1414bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)

  • branches/pap/ppStack/src/ppStackConvolve.c

    r25955 r26007  
    99#include "ppStack.h"
    1010#include "ppStackLoop.h"
     11
     12//#define TESTING
     13
    1114
    1215// Update the value of a concept
    … …  
    130133        ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config);
    131134        psFree(maskHeader);
    132         psImageCovarianceTransfer(readout->variance, readout->covariance);
    133135        ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config);
    134136#ifdef TESTING
    … …  
    139141            pmStackVisualPlotTestImage(readout->covariance->image, name);
    140142            psFree(name);
     143        }
     144        {
     145            int numCols = readout->image->numCols, numRows = readout->image->numRows;
     146            psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     147            for (int y = 0; y < numRows; y++) {
     148                for (int x = 0; x < numCols; x++) {
     149                    sn->data.F32[y][x] = readout->image->data.F32[y][x] /
     150                        sqrtf(readout->variance->data.F32[y][x]);
     151                }
     152            }
     153            psString name = NULL;
     154            psStringAppend(&name, "signoise_%d.fits", i);
     155            ppStackWriteImage(name, hdu->header, sn, config);
     156            psFree(name);
     157            psFree(sn);
    141158        }
    142159#endif
    … …  
    221238            numGood = 0;                    // Number of good images
    222239            for (int i = 0; i < num; i++) {
    223               if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
     240                if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    224241                    continue;
    225242                }
    … …  
    245262    // Correct chi^2 for renormalisation
    246263    psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms);
     264    for (int i = 0; i < num; i++) {
     265        psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n",
     266                 i, options->matchChi2->data.F32[i]);
     267    }
    247268    psFree(renorms);
    248269

  • branches/pap/ppStack/src/ppStackLoop.c

    r25964 r26007  
    130130        }
    131131        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
    132         if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, true, true)) {
     132        if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) {
    133133            psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination.");
    134134            psFree(stack);

  • branches/pap/ppStack/src/ppStackMatch.c

    r25959 r26007  
    1818#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
    1919
    20 #define TESTING                         // Enable debugging output
     20//#define TESTING                         // Enable debugging output
    2121
    2222#ifdef TESTING
    … …  
    167167    )
    168168{
     169#if 1
    169170    bool mdok; // Status of metadata lookups
    170171
    … …  
    192193    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
    193194
     195    psImageCovarianceTransfer(readout->variance, readout->covariance);
    194196    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
     197#else
     198    return true;
     199#endif
    195200}
    196201
    … …  
    212217    int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    213218
    214     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     219    psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
    215220    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    216221    psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
    … …  
    377382            }
    378383#endif
     384
     385            fprintf(stderr, "vf = %f\n", psImageCovarianceFactor(readout->covariance));
     386
    379387
    380388            if (threads > 0) {
    … …  
    516524            psFree(iter);
    517525            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     526            fprintf(stderr, "chi2 = %f ; vf = %f\n", sum/num, psImageCovarianceFactor(readout->covariance));
    518527        }
    519528
    … …  
    542551    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    543552    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    544       psWarning("Can't measure background for image.");
    545       psErrorClear();
     553        psWarning("Can't measure background for image.");
     554        psErrorClear();
    546555    } else {
    547       if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
    548         psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
    549                  psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
    550         (void)psBinaryOp(readout->image, readout->image, "-",
    551                          psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
    552       }
     556        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
     557            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
     558                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
     559            (void)psBinaryOp(readout->image, readout->image, "-",
     560                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
     561        }
    553562    }
    554563

  • branches/pap/ppStack/src/ppStackReadout.c

    r25968 r26007  
    116116
    117117    bool mdok;                          // Status of MD lookup
    118     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     118    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
    119119    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    120120    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
    … …  
    126126    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    127127
    128     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     128    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    129129    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     130    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     131    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    130132    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    131133    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    161163    }
    162164
    163     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
     165    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
    164166                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    165167        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    … …  
    217219    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    218220
    219     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     221    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    220222    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     223    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     224    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    221225    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    222226    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    267271    }
    268272
    269     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, iter, combineRej,
     273    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
    270274                        combineSys, combineDiscard, useVariance, safe, true)) {
    271275        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");

  • branches/pap/ppStack/src/ppStackReject.c

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING
     12//#define TESTING
    1313
    1414bool ppStackReject(ppStackOptions *options, pmConfig *config)
Note: See TracChangeset for help on using the changeset viewer.