IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5807


Ignore:
Timestamp:
Dec 17, 2005, 11:23:59 AM (21 years ago)
Author:
magnier
Message:
  • adding 3rd order polynomial for ApTrend
  • dropping old version of pmPSFtryMetric
  • renamed pmPSFtryMetric_Alt to pmPSFtryMetric
  • renamed try to psfTry
  • dropped psLibUtils.h, psModuleUtils.h
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c

    r5751 r5807  
    11# include <pslib.h>
    2 # include "psLibUtils.h"
    32# include "pmObjects.h"
    4 # include "psModulesUtils.h"
    53# include "pmPSF.h"
    64# include "pmPSFtry.h"
    … …  
    4139    type           = pmModelSetType (modelName);
    4240    test->psf      = pmPSFAlloc (type);
    43     test->sources  = psMemCopy(sources);
     41    test->sources  = psMemIncrRefCounter(sources);
    4442    test->modelFLT = psArrayAlloc (sources->n);
    4543    test->modelPSF = psArrayAlloc (sources->n);
    … …  
    7977    int Npsf = 0;
    8078
    81     pmPSFtry *try
    82     = pmPSFtryAlloc (sources, modelName);
     79    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
    8380
    8481    // stage 1:  fit an independent model (freeModel) to all sources
    8582    psTimerStart ("fit");
    86     for (int i = 0; i < try
    87                 ->sources->n; i++) {
    88 
    89             pmSource *source = try
    90                                    ->sources->data[i];
    91             pmModel  *model  = pmSourceModelGuess (source, try
    92                                                        ->psf->type);
    93             x = source->peak->x;
    94             y = source->peak->y;
    95 
    96             // set temporary object mask and fit object
    97             // fit model as FLT, not PSF
    98             psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    99             status = pmSourceFitModel (source, model, false);
    100             psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    101 
    102             // exclude the poor fits
    103             if (!status) {
    104                 try
    105                     ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
    106                 psFree (model);
    107                 continue;
    108             }
    109             try
    110                 ->modelFLT->data[i] = model;
    111             Nflt ++;
     83    for (int i = 0; i < psfTry->sources->n; i++) {
     84
     85        pmSource *source = psfTry->sources->data[i];
     86        pmModel  *model  = pmSourceModelGuess (source, psfTry->psf->type);
     87        x = source->peak->x;
     88        y = source->peak->y;
     89
     90        // set temporary object mask and fit object
     91        // fit model as FLT, not PSF
     92        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
     93        status = pmSourceFitModel (source, model, false);
     94        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     95
     96        // exclude the poor fits
     97        if (!status) {
     98            psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     99            psFree (model);
     100            continue;
    112101        }
     102        psfTry->modelFLT->data[i] = model;
     103        Nflt ++;
     104    }
    113105    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    114106    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    115107
    116108    // make this optional?
    117     // DumpModelFits (try->modelFLT, "modelsFLT.dat");
     109    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
    118110
    119111    // stage 2: construct a psf (pmPSF) from this collection of model fits
    120     pmPSFFromModels (try
    121                      ->psf, try
    122                          ->modelFLT, try
    123                              ->mask);
     112    pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
    124113
    125114    // stage 3: refit with fixed shape parameters
    126115    psTimerStart ("fit");
    127     for (int i = 0; i < try
    128                 ->sources->n; i++) {
    129             // masked for: bad model fit, outlier in parameters
    130             if (try
    131                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    132 
    133             pmSource *source = try
    134                                    ->sources->data[i];
    135             pmModel  *modelFLT = try
    136                                      ->modelFLT->data[i];
    137 
    138             // set shape for this model based on PSF
    139             pmModel *modelPSF = pmModelFromPSF (modelFLT, try
    140                                                     ->psf);
    141             x = source->peak->x;
    142             y = source->peak->y;
    143 
    144             psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    145             status = pmSourceFitModel (source, modelPSF, true);
    146 
    147             // skip poor fits
    148             if (!status) {
    149                 try
    150                     ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
    151                 psFree (modelPSF);
    152                 goto next_source;
    153             }
    154 
    155             // otherwise, save the resulting model
    156             try
    157                 ->modelPSF->data[i] = modelPSF;
    158 
    159             // XXX : use a different aperture radius from the fit radius?
    160             // XXX : use a different estimator for the local sky?
    161             // XXX : pass 'source' as input?
    162             if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
    163                 try
    164                     ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    165                 goto next_source;
    166             }
    167 
    168             try
    169                 ->metric->data.F64[i] = obsMag - fitMag;
    170             try
    171                 ->fitMag->data.F64[i] = fitMag;
    172             Npsf ++;
     116    for (int i = 0; i < psfTry->sources->n; i++) {
     117        // masked for: bad model fit, outlier in parameters
     118        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
     119            continue;
     120
     121        pmSource *source = psfTry->sources->data[i];
     122        pmModel  *modelFLT = psfTry->modelFLT->data[i];
     123
     124        // set shape for this model based on PSF
     125        pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
     126        x = source->peak->x;
     127        y = source->peak->y;
     128
     129        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
     130        status = pmSourceFitModel (source, modelPSF, true);
     131
     132        // skip poor fits
     133        if (!status) {
     134            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
     135            psFree (modelPSF);
     136            goto next_source;
     137        }
     138
     139        // otherwise, save the resulting model
     140        psfTry->modelPSF->data[i] = modelPSF;
     141
     142        // XXX : use a different aperture radius from the fit radius?
     143        // XXX : use a different estimator for the local sky?
     144        // XXX : pass 'source' as input?
     145        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
     146            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
     147            goto next_source;
     148        }
     149
     150        psfTry->metric->data.F64[i] = obsMag - fitMag;
     151        psfTry->fitMag->data.F64[i] = fitMag;
     152        Npsf ++;
    173153
    174154next_source:
    175             psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    176 
    177         }
     155        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     156
     157    }
    178158    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    179159    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    180160
    181161    // make this optional
    182     // DumpModelFits (try->modelPSF, "modelsPSF.dat");
     162    // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
    183163
    184164    // XXX this function wants aperture radius for pmSourcePhotometry
    185     pmPSFtryMetric_Alt (try
    186                         , RADIUS);
    187     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    188               modelName, try
    189                   ->psf->ApResid, try
    190                       ->psf->dApResid, try
    191                           ->psf->skyBias);
    192 
    193     return (try
    194            );
    195 }
    196 
    197 
    198 bool pmPSFtryMetric (pmPSFtry *try
    199                      , float RADIUS)
    200 {
    201 
    202     float dBin;
    203     int   nKeep, nSkip;
    204 
    205     // the measured (aperture - fit) magnitudes (dA == try->metric)
     165    pmPSFtryMetric (psfTry, RADIUS);
     166    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
     167              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
     168
     169    return (psfTry);
     170}
     171
     172bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
     173{
     174
     175    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    206176    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    207177    //     dA = dAo + dsky/flux
    … …  
    213183
    214184    // rflux = ten(0.4*fitMag);
    215     psVector *rflux = psVectorAlloc (try
    216                                      ->sources->n, PS_TYPE_F64);
    217     for (int i = 0; i < try
    218                 ->sources->n; i++) {
    219             if (try
    220                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    221             rflux->data.F64[i] = pow(10.0, 0.4*try
    222                                      ->fitMag->data.F64[i]);
    223         }
    224 
    225     // find min and max of (1/flux):
    226     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    227     psVectorStats (stats, rflux, NULL, try
    228                        ->mask, PSFTRY_MASK_ALL);
    229 
    230     // build binned versions of rflux, metric
    231     dBin = (stats->max - stats->min) / 10.0;
    232     psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    233     psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    234     psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    235     psFree (stats);
    236 
    237     psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    238 
    239     // group data in daBin bins, measure lower 50% mean
    240     for (int i = 0; i < daBin->n; i++) {
    241 
    242         psVector *tmp = psVectorAlloc (try
    243                                        ->sources->n, PS_TYPE_F64);
    244         tmp->n = 0;
    245 
    246         // accumulate data within bin range
    247         for (int j = 0; j < try
    248                     ->sources->n; j++) {
    249                 // masked for: bad model fit, outlier in parameters
    250                 if (try
    251                         ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
    252 
    253                 // skip points with extreme dA values
    254                 if (fabs(try
    255                          ->metric->data.F64[j]) > 0.5) continue;
    256 
    257                 // skip points outside of this bin
    258                 if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    259                     continue;
    260                 if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    261                     continue;
    262 
    263                 tmp->data.F64[tmp->n] = try
    264                                             ->metric->data.F64[j];
    265                 tmp->n ++;
    266             }
    267 
    268         // is this a valid point?
    269         maskB->data.U8[i] = 0;
    270         if (tmp->n < 2) {
    271             maskB->data.U8[i] = 1;
    272             psFree (tmp);
     185    psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     186    for (int i = 0; i < psfTry->sources->n; i++) {
     187        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
    273188            continue;
    274         }
    275 
    276         // dA values are contaminated with low outliers
    277         // measure statistics only on upper 50% of points
    278         // this would be easier if we could sort in reverse:
    279         //
    280         // psVectorSort (tmp, tmp);
    281         // tmp->n = 0.5*tmp->n;
    282         // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    283         // psVectorStats (stats, tmp, NULL, NULL, 0);
    284         // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    285 
    286         psVectorSort (tmp, tmp);
    287         nKeep = 0.5*tmp->n;
    288         nSkip = tmp->n - nKeep;
    289 
    290         psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    291         for (int j = 0; j < tmp2->n; j++) {
    292             tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    293         }
    294 
    295         stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    296         psVectorStats (stats, tmp2, NULL, NULL, 0);
    297         psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    298 
    299         daBin->data.F64[i] = stats->sampleMedian;
    300 
    301         psFree (stats);
    302         psFree (tmp);
    303         psFree (tmp2);
    304     }
    305 
    306     // linear fit to rfBin, daBin
    307     psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    308     psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    309     poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
    310 
    311     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    312     // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
    313 
    314     // XXX EAM : replace this when the above version is implemented
    315     // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
    316 
    317     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    318     psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    319 
    320     stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    321     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    322 
    323     try
    324         ->psf->ApResid = poly->coeff[0];
    325     try
    326         ->psf->dApResid = stats->clippedStdev;
    327     try
    328         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    329 
    330     psFree (rflux);
    331     psFree (rfBin);
    332     psFree (daBin);
    333     psFree (maskB);
    334     psFree (daBinFit);
    335     psFree (daResid);
    336     psFree (poly);
    337     psFree (stats);
    338     psFree (fitstat);
    339 
    340     return true;
    341 }
    342 
    343 bool pmPSFtryMetric_Alt (pmPSFtry *try
    344                          , float RADIUS)
    345 {
    346 
    347     // the measured (aperture - fit) magnitudes (dA == try->metric)
    348     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    349     //     dA = dAo + dsky/flux
    350     //   where flux is the flux of the star
    351     // we fit this trend to find the infinite flux aperture correction (dAo),
    352     //   the nominal sky bias (dsky), and the error on dAo
    353     // the values of dA are contaminated by stars with close neighbors in the aperture
    354     //   we use an outlier rejection to avoid this bias
    355 
    356     // rflux = ten(0.4*fitMag);
    357     psVector *rflux = psVectorAlloc (try
    358                                      ->sources->n, PS_TYPE_F64);
    359     for (int i = 0; i < try
    360                 ->sources->n; i++) {
    361             if (try
    362                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    363             rflux->data.F64[i] = pow(10.0, 0.4*try
    364                                      ->fitMag->data.F64[i]);
    365         }
    366 
    367     // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     189        rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
     190    }
     191
     192    // use 3hi/1lo sigma clipping on the rflux vs metric fit
    368193    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    369 
    370     // XXX EAM
    371194    stats->min = 1.0;
    372195    stats->max = 3.0;
    373196    stats->clipIter = 3;
    374197
    375     // linear clipped fit to rfBin, daBin
     198    // linear clipped fit of ApResid to rflux
    376199    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    377     poly = psVectorClipFitPolynomial1D (poly, stats, try
    378                                             ->mask, PSFTRY_MASK_ALL, try
    379                                                 ->metric, NULL, rflux);
    380     fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    381 
    382     try
    383         ->psf->ApResid = poly->coeff[0];
    384     try
    385         ->psf->dApResid = stats->sampleStdev;
    386     try
    387         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    388 
    389     fprintf (stderr, "*******************************************************************************\n");
    390 
    391     FILE *f;
    392     f = fopen ("apresid.dat", "w");
    393     if (f == NULL)
    394         psAbort ("pmPSFtry", "can't open output file");
    395 
    396     for (int i = 0; i < try
    397                 ->sources->n; i++) {
    398             fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
    399                          ->fitMag->data.F64[i], rflux->data.F64[i], try
    400                              ->metric->data.F64[i], try
    401                                  ->mask->data.U8[i]);
    402         }
    403     fclose (f);
     200    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, rflux);
     201    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     202
     203    psfTry->psf->ApTrend->coeff[0][0][0] = poly->coeff[0];
     204    psfTry->psf->ApTrend->coeff[0][0][1] = poly->coeff[1];
     205
     206    psfTry->psf->skyBias  = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     207    psfTry->psf->ApResid  = poly->coeff[0];
     208    psfTry->psf->dApResid = stats->sampleStdev;
    404209
    405210    psFree (rflux);
    … …  
    407212    psFree (stats);
    408213
    409     // psFree (daFit);
    410     // psFree (daResid);
    411 
    412214    return true;
    413215}
Note: See TracChangeset for help on using the changeset viewer.