IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15842


Ignore:
Timestamp:
Dec 14, 2007, 3:23:18 PM (18 years ago)
Author:
Paul Price
Message:

Making it so that the user can select the type of statistic to use when fitting a PSF.

Location:
trunk/psModules/src/objects
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r15705 r15842  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-11-28 00:43:27 $
     7 *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-12-15 01:23:18 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    221221
    222222    // use 3hi/3lo sigma clipping on the chisq fit
    223     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    224     stats->clipSigma = 3.0;
    225     stats->clipIter = 3;
     223    psStats *stats = options->stats;
    226224
    227225    // linear clipped fit of chisq trend vs flux
    228     bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, stats, mask, 0xff, chisq, NULL, flux);
    229     psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
    230 
    231     psFree (flux);
    232     psFree (stats);
    233     psFree (mask);
    234     psFree (chisq);
     226    bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, mask,
     227                                              0xff, chisq, NULL, flux);
     228    psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
     229    psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
     230
     231    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
     232              psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
     233
     234    psFree(flux);
     235    psFree(mask);
     236    psFree(chisq);
    235237
    236238    if (!result) {
    … …  
    247249
    248250    // XXX this function wants aperture radius for pmSourcePhotometry
    249     if (!pmPSFtryMetric (psfTry, options->radius)) {
     251    if (!pmPSFtryMetric (psfTry, options)) {
    250252        psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
    251253        return NULL;
    … …  
    258260}
    259261
    260 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
     262bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
    261263{
    262264    PS_ASSERT_PTR_NON_NULL(psfTry, false);
     265    PS_ASSERT_PTR_NON_NULL(options, false);
    263266    PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
     267
     268    float RADIUS = options->radius;
    264269
    265270    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    … …  
    286291    // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
    287292
    288     // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
    289     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    290     stats->clipSigma = 3.0;
    291     stats->clipIter = 3;
    292 
    293293    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    294294    // linear clipped fit of ApResid to r2rflux
    … …  
    297297
    298298    // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    299     bool result = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, psfTry->metricErr, r2rflux);
     299    bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
     300                                              psfTry->metric, psfTry->metricErr, r2rflux);
    300301    if (!result) {
    301302        psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");
    … …  
    303304        psFree(poly);
    304305        psFree(r2rflux);
    305         psFree(stats);
    306306
    307307        return false;
    308308    }
    309     psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; keeping %ld of %ld psf stars\n", poly->coeff[0], stats->robustStdev, stats->clippedNvalues, psfTry->sources->n);
     309    psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
     310    psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
     311              psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
     312
     313
    310314
    311315    // XXX test dump of fitted model (dump when tracing?)
    … …  
    335339    psfTry->psf->skyBias  = poly->coeff[1];
    336340    psfTry->psf->ApResid  = poly->coeff[0];
    337     psfTry->psf->dApResid = stats->robustStdev;
     341    psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);
    338342
    339343    psFree (r2rflux);
    340344    psFree (poly);
    341     psFree (stats);
    342345
    343346    return true;
    … …  
    720723    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    721724    int nGroup = PS_MAX (3*Nx*Ny, 10);
    722     pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup);
     725    pmPSFShapeParamsErrors (errorFloor, mag, e0res, e1res, e2res, mask, nGroup,
     726                            psStatsStdevOption(psf->psfTrendStats->options));
    723727
    724728    *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
    … …  
    739743
    740744// calculate the minimum scatter of the parameters
    741 bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup) {
    742 
    743     psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     745bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
     746                            psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
     747{
     748
     749    psStats *statsS = psStatsAlloc(stdevOpt);
    744750
    745751    // measure the trend in bins with 10 values each; if < 10 total, use them all
    … …  
    778784        psStatsInit (statsS);
    779785        psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff);
    780         dEsquare += PS_SQR(statsS->robustStdev);
     786        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    781787
    782788        psStatsInit (statsS);
    783789        psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff);
    784         dEsquare += PS_SQR(statsS->robustStdev);
     790        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    785791
    786792        psStatsInit (statsS);
    787793        psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff);
    788         dEsquare += PS_SQR(statsS->robustStdev);
     794        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    789795
    790796        dSo->data.F32[i] = sqrt(dEsquare);
    … …  
    795801    psFree (mkSubset);
    796802
    797     psStats *stats = psStatsAlloc (PS_STAT_MIN);
    798     psVectorStats (stats, dSo, NULL, NULL, 0);
    799 
    800     *errorFloor = stats->min;
    801 
    802     psFree (stats);
    803     psFree (index);
    804 
     803    psStats *minStats = psStatsAlloc(PS_STAT_MIN);
     804    psVectorStats(minStats, dSo, NULL, NULL, 0);
     805    *errorFloor = minStats->min;
     806    psFree(minStats);
     807
     808    psFree(index);
    805809    psFree (dSo);
    806810

  • trunk/psModules/src/objects/pmPSFtry.h

    r15697 r15842  
    66 * @author EAM, IfA
    77 *
    8  * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2007-11-27 03:14:57 $
     8 * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2007-12-15 01:23:18 $
    1010 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 */
    … …  
    9999bool pmPSFtryMetric(
    100100    pmPSFtry *psfTry,                  ///< Add comment.
    101     float RADIUS                       ///< Add comment.
     101    pmPSFOptions *options              ///< PSF fitting options
    102102);
    103103
    … …  
    129129bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
    130130bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *errorFloor, float *errorTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz);
    131 bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup);
     131bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt);
    132132
    133133/// @}

  • trunk/psModules/src/objects/pmTrend2D.h

    r15696 r15842  
    55 * @author EAM, IfA
    66 *
    7  * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2007-11-27 03:08:38 $
     7 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2007-12-15 01:23:18 $
    99 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1010 */
    … …  
    1212# ifndef PM_TREND_2D_H
    1313# define PM_TREND_2D_H
     14
     15#include <pslib.h>
    1416
    1517/// @addtogroup Objects Object Detection / Analysis Functions
    … …  
    2628    psPolynomial2D *poly;
    2729    psImageMap *map;
    28     psStats *stats;
     30    psStats *stats;                     // Statistics for clipped fitting
     31    psStatsOptions singleMean, singleStdev; // Staistics for mean and stdev when there's a single pixel
    2932    pmTrend2DMode mode; // POLY_ORD, POLY_CHEB, MAP
    3033} pmTrend2D;
Note: See TracChangeset for help on using the changeset viewer.