IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35455


Ignore:
Timestamp:
Apr 30, 2013, 6:10:41 PM (13 years ago)
Author:
watersc1
Message:

Changes to stacking code:

  • New input FWHM clipping algorithm using mixture models
  • New recipes/reductions that further reduce FWHM clipping
  • pmSubtractionRejectStamps now fits flux/chi2 values in log-log space to prevent outliers biasing the fit (and causing images to be rejected due to a small number of bad stamps)
  • pmPSF now correctly chooses between useReff forms as appropriate
  • pmPSF change fixes circularization issues in pmReadoutFakeFromSources
  • pmPSFEnvelope somewhat more clear about circularization
  • Changed debugging in pmSubtractionEquation to add information.

I do not see any conflicts with the ppSub code.

Location:
trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippconfig/gpc1/ppStack.config

    r35399 r35455  
    6363    PSF.INPUT.THRESH        F32   5.0
    6464    PSF.INPUT.ASYMMETRY     F32   0.2
     65   MATCH.REJ        F32    4.0
     66   SAFE          BOOL   F
     67
     68END
     69
     70STACK_THREEPI_ALT  METADATA
     71    OUTPUT.NOCOMP           BOOL  TRUE  ## MD test mod
     72    OUTPUT.LOGFLUX          BOOL  FALSE ## MD test mod
     73    STACK.TYPE              STR   DEEP_STACK ## MD test mod
     74    PSF.INPUT.MAX           F32   10.0  # input cut more restrictive now 12->10 pix, limit target PSF in convolved stacks
     75    PSF.INPUT.CLIP.NSIGMA   F32   1.5   # additional cut if majority of inputs much smaller than INPUT.MAX
     76    PSF.INPUT.THRESH        F32   10.0
     77    PSF.INPUT.ASYMMETRY     F32   0.2
     78   MATCH.REJ        F32    4.0
     79   SAFE          BOOL   F
     80
    6581END
    6682

  • trunk/ippconfig/recipes/ppStack.config

    r35394 r35455  
    160160STACK_THREEPI  METADATA
    161161END
     162STACK_THREEPI_ALT  METADATA
     163END
    162164
    163165STACK_ALLDEEP   METADATA

  • trunk/ippconfig/recipes/reductionClasses.mdc

    r35337 r35455  
    538538END
    539539
     540THREEPI_STACK_ALT     METADATA
     541      STACK_PPSTACK   STR      STACK_THREEPI_ALT
     542      STACK_PPSUB     STR      STACK
     543      STACK_PSPHOT    STR      STACK
     544END
     545
     546
    540547# quick stacks
    541548QUICKSTACK            METADATA
    … …  
    632639PRESERVE_BG     METADATA
    633640        CHIP_PPIMAGE    STR     CHIP_PRESERVE_BG
     641        CHIP_PSPHOT     STR     CHIP
     642        WARP_PSWARP     STR     WARP
     643        STACK_PPSTACK   STR     STACK
     644        STACK_PPSUB     STR     STACK
     645        STACK_PSPHOT    STR     STACK
     646        DIFF_PPSUB      STR     DIFF
     647        DIFF_PSPHOT     STR     DIFF
     648        JPEG_BIN1       STR     PPIMAGE_J1
     649        JPEG_BIN2       STR     PPIMAGE_J2
     650        FAKEPHOT        STR     FAKEPHOT
     651        ADDSTAR         STR     ADDSTAR
     652        PSASTRO         STR     DEFAULT_RECIPE
     653        BACKGROUND_PPBACKGROUND STR     BACKGROUND
     654        BACKGROUND_PSWARP       STR     BACKGROUND
     655END
     656
     657# Background removal, but then add it back in
     658RESTORE_BG      METADATA
     659        CHIP_PPIMAGE    STR     CHIP_RESTORE_BG
    634660        CHIP_PSPHOT     STR     CHIP
    635661        WARP_PSWARP     STR     WARP

  • trunk/ppStack/src/ppStackConvolve.c

    r34372 r35455  
    403403                    continue;
    404404                }
    405                 if ((options->matchChi2->data.F32[i] > thresh) || ! isfinite(options->matchChi2->data.F32[i])) {
     405                if ((options->matchChi2->data.F32[i] > thresh) ||
     406                    ! isfinite(options->matchChi2->data.F32[i])) {
    406407                    numRej++;
    407408                    options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;

  • trunk/ppStack/src/ppStackMatch.c

    r35383 r35455  
    467467                sum += kernels->rms;
    468468                num++;
     469                //              psLogMsg("ppStack", PS_LOG_INFO, "KERNMATCHDUMP %d %d %g %g %g\n",
     470                //                       index,num,kernels->mean,kernels->rms,sum);
    469471            }
    470472            psFree(iter);
    471473            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     474
    472475        }
    473476

  • trunk/ppStack/src/ppStackPrepare.c

    r35400 r35455  
    334334    //    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
    335335    //    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
    336 
    337336    //    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
    338    
    339     //    unwrittenGMMfunction(options->inputSeeing, options->inputMask, 0xff,
    340     //           &Punimodal,
    341     //           &m1,&s1,&pi1,
    342     //           &m2,&s2,&pi2);
    343337
    344338    // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
    … …  
    358352    if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
    359353    if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
    360 
     354    psLogMsg("ppStack",PS_LOG_INFO,
     355             "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
     356             limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
     357             means->data.F32[1],S->data.F32[1],pi->data.F32[1],
     358             Punimodal);
    361359    for (int i = 0; i < num; i++) {
    362360      if (options->inputSeeing->data.F32[i] > limit) {
    … …  
    404402    if (options->convolve) {
    405403        options->psf = ppStackPSF(config, numCols, numRows, psfs, options);
    406         psFree(psfs);
     404        psFree(psfs);
    407405        if (!options->psf) {
    408406#if 1

  • trunk/psModules/src/camera/pmReadoutFake.c

    r34403 r35455  
    152152            return false;
    153153        }
    154 
     154       
    155155        psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
    156156                fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r34403 r35455  
    3737
    3838
    39 //#define TESTING                         // Enable test output
     39// #define TESTING                         // Enable test output
    4040// #define PEAK_NORM                       // Normalise peaks?
    4141#define PEAK_FLUX 1.0e4                 // Peak flux for each source
    … …  
    190190        psf->residuals = NULL;
    191191        pmModelClassSetLimits(PM_MODEL_LIMITS_MODERATE);
     192        psLogMsg("psModules",PS_LOG_INFO,"Matching Input %d",i);
     193#define CIRCULARIZE true
    192194        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
    193                                       NAN, radius, true, false)) {
     195                                      NAN, radius, CIRCULARIZE, false)) {
    194196            psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
    195197            psFree(envelope);

  • trunk/psModules/src/imcombine/pmStack.c

    r35165 r35455  
    4141/* #define TEST_X 5745                       // x coordinate to examine */
    4242/* #define TEST_Y 5331                       // y coordinate to examine */
    43 #define TEST_X 25
    44 #define TEST_Y 25
     43// #define TEST_X 972
     44// #define TEST_Y 3213
     45#define TEST_X 3289
     46#define TEST_Y 4810
    4547#define TEST_RADIUS 0.5                 // Radius to examine
    4648# endif
    … …  
    752754    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    753755        for (int i = 0; i < numGood; i++) {
    754             fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
    755                     i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    756                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
    757                     pixelSuspects->data.U8[i], badMaskBits, suspectMaskBits, *badMask, *goodMask);
     756            fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%g) %g %f %d %x %x -> %x %x\n",
     757                    i, x, y, pixelSources->data.U16[i],
     758                    pixelData->data.F32[i], pixelVariances->data.F32[i],
     759                    addVariance->data.F32[i],
     760                    pixelWeights->data.F32[i], pixelExps->data.F32[i],
     761                    pixelSuspects->data.U8[i],
     762                    badMaskBits, suspectMaskBits,
     763                    *badMask, *goodMask);
    758764        }
    759765    }

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r33425 r35455  
    3737#define USE_KERNEL_ERR                  // Use kernel error image?
    3838#define NUM_COVAR_POS 5                 // Number of positions for covariance calculation
     39#define USE_LOGFIT_REJECT
    3940
    4041// XXX we need to pass these fwhm values elsewhere.  These should go on one of the structure, but
    … …  
    966967    psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
    967968
     969#ifdef USE_LOGFIT_REJECT
     970    // CZW: Since flux and chisq can range over many orders of magnitude, use a log-log fit to calculate the
     971    // RMS scatter.  This should be more robust against outliers that can impact the fit quality.
     972    psVector *logchisq = psVectorAlloc(match->chisq->n,PS_TYPE_F32);
     973    psVector *logflux  = psVectorAlloc(match->fluxes->n,PS_TYPE_F32);
     974    for (int qq = 0; qq < match->chisq->n; qq++) {
     975      if ((match->chisq->data.F32[qq] > 0)&&
     976          (match->fluxes->data.F32[qq] > 0)) {
     977        logchisq->data.F32[qq] = log10(match->chisq->data.F32[qq]);
     978        logflux->data.F32[qq]  = log10(match->fluxes->data.F32[qq]);
     979      }
     980      else {
     981        // Ignore negative values.  Maybe we could do an offset or something, but not today.
     982        match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] |= PM_SUBTRACTION_STAMP_REJECTED;
     983      }
     984      //      if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) {
     985      //        psLogMsg("psModules.imcombine", PS_LOG_INFO, "RRRRRR: %d %g %g %g %g\n",
     986      //                 qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq],
     987      //                 0.0,0.0);
     988      // }
     989
     990    }
     991    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, logchisq, NULL, logflux);
     992    if (!result) {
     993        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
     994        psFree(model);
     995        psFree(stats);
     996        return -1;
     997    }
     998
     999    // However, since we've done a log-log fit, we can't rely on the statistics in the stats object
     1000    // To accurately represent the distribution.  I've tried to massage them (as dlog10(X) = abs(1/(X log(10))) dX),
     1001    // but ended up deciding to just manually calculate the residual, and then pass that to a robust stats object.
     1002    psStats *residStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     1003    psVector *residData = psVectorAlloc(match->chisq->n,PS_TYPE_F32);
     1004    int counter = 0;
     1005    for (int qq = 0; qq < match->chisq->n; qq++) {
     1006      if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) {
     1007        double vv = pow(10,
     1008                        model->coeff[0] +
     1009                        model->coeff[1] * log10(match->fluxes->data.F32[qq]) +
     1010                        model->coeff[2] * pow(log10(match->fluxes->data.F32[qq]),2));
     1011        residData->data.F32[qq] = match->chisq->data.F32[qq] - vv;
     1012        counter++;
     1013        //                      psLogMsg("psModules.imcombine", PS_LOG_INFO, "SSSSS: %d %g %g %g %g\n",
     1014        //                       qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq],
     1015        //                       vv,residData->data.F32[qq]);
     1016      }
     1017    }
     1018    psVectorStats(residStats,residData,NULL,match->stampMask,0xff);
     1019    if (isnan(residStats->robustMedian)) {
     1020        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
     1021        psFree(model);
     1022        psFree(stats);
     1023        psFree(logchisq);
     1024        psFree(logflux);
     1025        psFree(residStats);
     1026        psFree(residData);
     1027
     1028        return -1;
     1029    }
     1030
     1031    kernels->mean = residStats->robustMedian;
     1032    kernels->rms =  residStats->robustStdev;
     1033    kernels->numStamps = counter;
     1034
     1035    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", residStats->robustMedian,residStats->robustStdev);
     1036    psFree(logchisq);
     1037    psFree(logflux);
     1038    psFree(residStats);
     1039    psFree(residData);
     1040
     1041#else
     1042    // CZW: Otherwise, use the original linear fit code.
    9681043    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes);
    9691044    if (!result) {
    … …  
    9791054        return -1;
    9801055    }
    981 
    9821056    kernels->mean = stats->sampleMean;
    983     kernels->rms = stats->sampleStdev;
     1057    kernels->rms =  stats->sampleStdev;
    9841058    kernels->numStamps = stats->clippedNvalues;
     1059    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
     1060#endif
    9851061
    9861062    psLogMsg ("pmPSFtry", 4, "chisq vs flux model: %e + %e flux + %e flux^2\n", model->coeff[0], model->coeff[1], model->coeff[2]);
    987     psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
    9881063    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",  kernels->numStamps, kernels->mean, kernels->rms);
    9891064

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r33089 r35455  
    751751    stamps->normValue2 = stats->robustMedian;
    752752
    753     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
     753    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f) %ld\n", stamps->normValue, stamps->normValue2,norms->n);
    754754
    755755    psFree(stats);
    … …  
    790790        }
    791791    }
    792 
     792    //    psLogMsg("psModules.imcombine",PS_LOG_INFO, "stampNorm: %f %f %d %d %d %f %f",
     793    //       stamp->x,stamp->y,footprint,normWindow1,normWindow2,
     794    //       normI1,normI2);
    793795    stamp->norm = normI2 / normI1;
    794796    stamp->normI1 = normI1;

  • trunk/psModules/src/objects/pmPSF.c

    r34403 r35455  
    342342    PS_ASSERT_PTR_NON_NULL(modelPar, axes);
    343343
    344     bool useReff = true;
     344    bool useReff = false;
    345345    useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
    346346    useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
    … …  
    350350        shape.sx  = modelPar[PM_PAR_SXX];
    351351        shape.sy  = modelPar[PM_PAR_SYY];
    352         shape.sxy = modelPar[PM_PAR_SXY];
     352        shape.sxy = modelPar[PM_PAR_SXY]; // This is potentially wrong?
    353353    } else {
    354354        shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
    … …  
    384384    psEllipseShape shape = psEllipseAxesToShape (axes);
    385385
    386     bool useReff = true;
    387     useReff |= pmModelClassGetType ("PS_MODEL_SERSIC");
    388     useReff |= pmModelClassGetType ("PS_MODEL_DEV");
    389     useReff |= pmModelClassGetType ("PS_MODEL_EXP");
     386    bool useReff = false;
     387    useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC"));
     388    useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV"));
     389    useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP"));
    390390
    391391    if (useReff) {
    392392        modelPar[PM_PAR_SXX] = shape.sx;
    393393        modelPar[PM_PAR_SYY] = shape.sy;
    394         modelPar[PM_PAR_SXY] = shape.sxy;
     394        modelPar[PM_PAR_SXY] = shape.sxy; // This is potentially wrong?
    395395    } else {
    396396        modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
Note: See TracChangeset for help on using the changeset viewer.