IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 17, 2006, 8:10:08 AM (20 years ago)
Author:
magnier
Message:

fixed up conflicts

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6511 r6873  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-03-04 01:01:33 $
     7 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-04-17 18:10:08 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1313
    1414# include <pslib.h>
    15 # include "pmObjects.h"
    16 # include "pmPSF.h"
    17 # include "pmPSFtry.h"
    18 # include "pmModelGroup.h"
     15#include "pmHDU.h"
     16#include "pmFPA.h"
     17#include "pmPeaks.h"
     18#include "pmMoments.h"
     19#include "pmModel.h"
     20#include "pmSource.h"
     21#include "pmGrowthCurve.h"
     22#include "pmPSF.h"
     23#include "pmPSFtry.h"
     24#include "pmModelGroup.h"
     25#include "pmSourceFitModel.h"
     26#include "pmSourcePhotometry.h"
    1927
    2028// ********  pmPSFtry functions  **************************************************
    … …  
    3341    psFree (test->psf);
    3442    psFree (test->sources);
    35     psFree (test->modelFLT);
     43    psFree (test->modelEXT);
    3644    psFree (test->modelPSF);
    3745    psFree (test->metric);
    … …  
    4250
    4351// allocate a pmPSFtry based on the desired sources and the model (identified by name)
    44 pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName)
     52pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors)
    4553{
    4654
    … …  
    5159    // XXX probably need to increment ref counter
    5260    type           = pmModelSetType (modelName);
    53     test->psf      = pmPSFAlloc (type);
     61    test->psf      = pmPSFAlloc (type, poissonErrors);
    5462    test->sources  = psMemIncrRefCounter(sources);
    55     test->modelFLT = psArrayAlloc (sources->n);
     63    test->modelEXT = psArrayAlloc (sources->n);
    5664    test->modelPSF = psArrayAlloc (sources->n);
    5765    test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
    … …  
    5967    test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
    6068
    61     test->modelFLT->n = test->modelFLT->nalloc;
    62     test->modelPSF->n = test->modelPSF->nalloc;
    63     test->metric->n   = test->metric->nalloc;
    64     test->fitMag->n   = test->fitMag->nalloc;
    65     test->mask->n     = test->mask->nalloc;
    66 
    67     for (int i = 0; i < test->modelFLT->n; i++) {
     69    for (int i = 0; i < test->modelEXT->n; i++) {
    6870        test->mask->data.U8[i]  = 0;
    69         test->modelFLT->data[i] = NULL;
     71        test->modelEXT->data[i] = NULL;
    7072        test->modelPSF->data[i] = NULL;
    7173        test->metric->data.F64[i] = 0;
    … …  
    8688// mask values indicate the reason the source was rejected:
    8789
    88 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
     90pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors)
    8991{
    9092    bool status;
    … …  
    9395    float x;
    9496    float y;
    95     int Nflt = 0;
     97    int Next = 0;
    9698    int Npsf = 0;
    9799
    98     pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
     100    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors);
    99101
    100102    // stage 1:  fit an independent model (freeModel) to all sources
    … …  
    108110
    109111        // set temporary object mask and fit object
    110         // fit model as FLT, not PSF
    111         psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    112         status = pmSourceFitModel (source, model, false);
    113         psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     112        // fit model as EXT, not PSF
     113
     114        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
     115        status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT);
     116        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED);
    114117
    115118        // exclude the poor fits
    116119        if (!status) {
    117             psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     120            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
    118121            psFree (model);
    119122            continue;
    120123        }
    121         psfTry->modelFLT->data[i] = model;
    122         Nflt ++;
    123     }
    124     psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    125     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    126 
    127     // make this optional?
    128     // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
     124        psfTry->modelEXT->data[i] = model;
     125        Next ++;
     126    }
     127    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     128    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n);
    129129
    130130    // stage 2: construct a psf (pmPSF) from this collection of model fits
    131     pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
     131    pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);
    132132
    133133    // stage 3: refit with fixed shape parameters
    … …  
    139139
    140140        pmSource *source = psfTry->sources->data[i];
    141         pmModel  *modelFLT = psfTry->modelFLT->data[i];
     141        pmModel  *modelEXT = psfTry->modelEXT->data[i];
    142142
    143143        // set shape for this model based on PSF
    144         pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
     144        pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
    145145        x = source->peak->x;
    146146        y = source->peak->y;
    147147
    148         psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    149         status = pmSourceFitModel (source, modelPSF, true);
     148        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
     149        status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF);
    150150
    151151        // skip poor fits
    … …  
    162162        // XXX : use a different estimator for the local sky?
    163163        // XXX : pass 'source' as input?
    164         if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
     164        if (!pmSourcePhotometryModel (&fitMag, modelPSF)) {
     165            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
     166            goto next_source;
     167        }
     168        if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) {
    165169            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    166170            goto next_source;
    … …  
    172176
    173177next_source:
    174         psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    175 
    176     }
     178        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED);
     179
     180    }
     181    psfTry->psf->nPSFstars = Npsf;
     182
    177183    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    178184    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    179185
    180     // make this optional
    181     // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
     186    // measure the chi-square trend as a function of flux (PAR[1])
     187    // this should be linear for Poisson errors and quadratic for constant sky errors
     188    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     189    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     190    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK);
     191
     192    // write sources with models first
     193    for (int i = 0; i < psfTry->sources->n; i++) {
     194        pmModel *model = psfTry->modelPSF->data[i];
     195        if (model == NULL) {
     196            flux->data.F64[i] = 0.0;
     197            chisq->data.F64[i] = 0.0;
     198            mask->data.U8[i] = 1;
     199        } else {
     200            flux->data.F64[i] = model->params->data.F32[1];
     201            chisq->data.F64[i] = model[0].chisq/model[0].nDOF;
     202            mask->data.U8[i] = 0;
     203        }
     204    }
     205    // use 3hi/3lo sigma clipping on the r2rflux vs metric fit
     206    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     207    stats->clipSigma = 3.0;
     208    stats->clipIter = 3;
     209    // linear clipped fit of ApResid to r2rflux
     210    psVectorClipFitPolynomial1D (psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
     211    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
     212        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), psfTry->psf->ChiTrend->coeffErr[i]);
     213    }
     214    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     215    psFree (stats);
     216    psFree (flux);
     217    psFree (chisq);
    182218
    183219    // XXX this function wants aperture radius for pmSourcePhotometry
    184     pmPSFtryMetric_Alt (psfTry, RADIUS);
    185     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    186               modelName,
    187               psfTry->psf->ApResid,
    188               psfTry->psf->dApResid,
    189               psfTry->psf->skyBias);
     220    pmPSFtryMetric (psfTry, RADIUS);
     221    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
     222              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    190223
    191224    return (psfTry);
    192225}
    193226
    194 
    195227bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
    196228{
    197 
    198     float dBin;
    199     int   nKeep, nSkip;
    200 
    201229    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    202230    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    … …  
    208236    //   we use an outlier rejection to avoid this bias
    209237
    210     // rflux = ten(0.4*fitMag);
    211     psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    212     rflux->n = rflux->nalloc;
     238    // r2rflux = radius^2 * ten(0.4*fitMag);
     239    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    213240    for (int i = 0; i < psfTry->sources->n; i++) {
    214241        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
    215242            continue;
    216         rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
    217     }
    218 
    219     // find min and max of (1/flux):
    220     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    221     psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
    222 
    223     // build binned versions of rflux, metric
    224     dBin = (stats->max - stats->min) / 10.0;
    225     psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    226     psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    227     psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    228     daBin->n = daBin->nalloc;
    229     maskB->n = maskB->nalloc;
    230     psFree (stats);
    231 
    232     psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    233 
    234     // group data in daBin bins, measure lower 50% mean
    235     for (int i = 0; i < daBin->n; i++) {
    236 
    237         psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    238         tmp->n = 0;
    239 
    240         // accumulate data within bin range
    241         for (int j = 0; j < psfTry->sources->n; j++) {
    242             // masked for: bad model fit, outlier in parameters
    243             if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
    244                 continue;
    245 
    246             // skip points with extreme dA values
    247             if (fabs(psfTry->metric->data.F64[j]) > 0.5)
    248                 continue;
    249 
    250             // skip points outside of this bin
    251             if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    252                 continue;
    253             if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    254                 continue;
    255 
    256             tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
    257             tmp->n ++;
    258         }
    259 
    260         // is this a valid point?
    261         maskB->data.U8[i] = 0;
    262         if (tmp->n < 2) {
    263             maskB->data.U8[i] = 1;
    264             psFree (tmp);
    265             continue;
    266         }
    267 
    268         // dA values are contaminated with low outliers
    269         // measure statistics only on upper 50% of points
    270         // this would be easier if we could sort in reverse:
    271         //
    272         // psVectorSort (tmp, tmp);
    273         // tmp->n = 0.5*tmp->n;
    274         // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    275         // psVectorStats (stats, tmp, NULL, NULL, 0);
    276         // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    277 
    278         psVectorSort (tmp, tmp);
    279         nKeep = 0.5*tmp->n;
    280         nSkip = tmp->n - nKeep;
    281 
    282         psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    283         tmp2->n = tmp2->nalloc;
    284         for (int j = 0; j < tmp2->n; j++) {
    285             tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    286         }
    287 
    288         stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    289         psVectorStats (stats, tmp2, NULL, NULL, 0);
    290         psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    291 
    292         daBin->data.F64[i] = stats->sampleMedian;
    293 
    294         psFree (stats);
    295         psFree (tmp);
    296         psFree (tmp2);
    297     }
    298 
    299     // linear fit to rfBin, daBin
    300     psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    301     psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    302     poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
    303 
    304     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    305     // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
    306 
    307     // XXX EAM : replace this when the above version is implemented
    308     // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
    309 
    310     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    311     psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    312 
    313     stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    314     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    315 
    316     psfTry->psf->ApResid = poly->coeff[0];
    317     psfTry->psf->dApResid = stats->clippedStdev;
    318     psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    319 
    320     psFree (rflux);
    321     psFree (rfBin);
    322     psFree (daBin);
    323     psFree (maskB);
    324     psFree (daBinFit);
    325     psFree (daResid);
    326     psFree (poly);
    327     psFree (stats);
    328     psFree (fitstat);
    329 
    330     return true;
    331 }
    332 
    333 bool pmPSFtryMetric_Alt (pmPSFtry *try
    334                          , float RADIUS)
    335 {
    336 
    337     // the measured (aperture - fit) magnitudes (dA == try->metric)
    338     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    339     //     dA = dAo + dsky/flux
    340     //   where flux is the flux of the star
    341     // we fit this trend to find the infinite flux aperture correction (dAo),
    342     //   the nominal sky bias (dsky), and the error on dAo
    343     // the values of dA are contaminated by stars with close neighbors in the aperture
    344     //   we use an outlier rejection to avoid this bias
    345 
    346     // rflux = ten(0.4*fitMag);
    347     psVector *rflux = psVectorAlloc (try
    348                                      ->sources->n, PS_TYPE_F64);
    349     rflux->n = rflux->nalloc;
    350     for (int i = 0; i < try
    351                 ->sources->n; i++) {
    352             if (try
    353                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    354             rflux->data.F64[i] = pow(10.0, 0.4*try
    355                                      ->fitMag->data.F64[i]);
    356         }
    357 
    358     // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     243        r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
     244    }
     245
     246    // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
    359247    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    360 
    361     // XXX EAM
    362248    stats->min = 1.0;
    363249    stats->max = 3.0;
    364250    stats->clipIter = 3;
    365251
    366     // linear clipped fit to rfBin, daBin
    367     psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    368     poly = psVectorClipFitPolynomial1D (poly, stats, try
    369                                             ->mask, PSFTRY_MASK_ALL, try
    370                                                 ->metric, NULL, rflux);
    371     fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    372 
    373     try
    374         ->psf->ApResid = poly->coeff[0];
    375     try
    376         ->psf->dApResid = stats->sampleStdev;
    377     try
    378         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    379 
    380     fprintf (stderr, "*******************************************************************************\n");
    381 
    382     FILE *f;
    383     f = fopen ("apresid.dat", "w");
    384     if (f == NULL)
    385         psAbort ("pmPSFtry", "can't open output file");
    386 
    387     for (int i = 0; i < try
    388                 ->sources->n; i++) {
    389             fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
    390                          ->fitMag->data.F64[i], rflux->data.F64[i], try
    391                              ->metric->data.F64[i], try
    392                                  ->mask->data.U8[i]);
    393         }
    394     fclose (f);
    395 
    396     psFree (rflux);
     252    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
     253    // linear clipped fit of ApResid to r2rflux
     254    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     255    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux);
     256    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     257
     258    pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS);
     259    psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
     260    psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;
     261
     262    psfTry->psf->skyBias  = poly->coeff[1];
     263    psfTry->psf->ApResid  = poly->coeff[0];
     264    psfTry->psf->dApResid = stats->sampleStdev;
     265
     266    psFree (r2rflux);
    397267    psFree (poly);
    398268    psFree (stats);
    399269
    400     // psFree (daFit);
    401     // psFree (daResid);
    402 
    403270    return true;
    404271}
     272
     273/*
     274  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
     275  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
     276  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
     277*/
     278
Note: See TracChangeset for help on using the changeset viewer.