IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 17, 2006, 7:13:42 AM (20 years ago)
Author:
magnier
Message:

bulk merge of eam_rel9_p0 onto this branch

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5844 r6448  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-12-24 01:24:32 $
     7 *  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-02-17 17:13:42 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3333    psFree (test->psf);
    3434    psFree (test->sources);
    35     psFree (test->modelFLT);
     35    psFree (test->modelEXT);
    3636    psFree (test->modelPSF);
    3737    psFree (test->metric);
    … …  
    5353    test->psf      = pmPSFAlloc (type);
    5454    test->sources  = psMemIncrRefCounter(sources);
    55     test->modelFLT = psArrayAlloc (sources->n);
     55    test->modelEXT = psArrayAlloc (sources->n);
    5656    test->modelPSF = psArrayAlloc (sources->n);
    5757    test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
    … …  
    5959    test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
    6060
    61     for (int i = 0; i < test->modelFLT->n; i++) {
     61    for (int i = 0; i < test->modelEXT->n; i++) {
    6262        test->mask->data.U8[i]  = 0;
    63         test->modelFLT->data[i] = NULL;
     63        test->modelEXT->data[i] = NULL;
    6464        test->modelPSF->data[i] = NULL;
    6565        test->metric->data.F64[i] = 0;
    … …  
    8787    float x;
    8888    float y;
    89     int Nflt = 0;
     89    int Next = 0;
    9090    int Npsf = 0;
    9191
    … …  
    102102
    103103        // set temporary object mask and fit object
    104         // fit model as FLT, not PSF
     104        // fit model as EXT, not PSF
     105
    105106        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    106107        status = pmSourceFitModel (source, model, false);
    … …  
    109110        // exclude the poor fits
    110111        if (!status) {
    111             psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     112            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
    112113            psFree (model);
    113114            continue;
    114115        }
    115         psfTry->modelFLT->data[i] = model;
    116         Nflt ++;
    117     }
    118     psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    119     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    120 
    121     // make this optional?
    122     // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
     116        psfTry->modelEXT->data[i] = model;
     117        Next ++;
     118    }
     119    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     120    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n);
    123121
    124122    // stage 2: construct a psf (pmPSF) from this collection of model fits
    125     pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
     123    pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);
    126124
    127125    // stage 3: refit with fixed shape parameters
    … …  
    133131
    134132        pmSource *source = psfTry->sources->data[i];
    135         pmModel  *modelFLT = psfTry->modelFLT->data[i];
     133        pmModel  *modelEXT = psfTry->modelEXT->data[i];
    136134
    137135        // set shape for this model based on PSF
    138         pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
     136        pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
    139137        x = source->peak->x;
    140138        y = source->peak->y;
    … …  
    169167
    170168    }
     169    psfTry->psf->nPSFstars = Npsf;
     170
    171171    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    172172    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    … …  
    176176
    177177    // XXX this function wants aperture radius for pmSourcePhotometry
    178     pmPSFtryMetric_Alt (psfTry, RADIUS);
    179     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    180               modelName,
    181               psfTry->psf->ApResid,
    182               psfTry->psf->dApResid,
    183               psfTry->psf->skyBias);
     178    pmPSFtryMetric (psfTry, RADIUS);
     179    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
     180              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    184181
    185182    return (psfTry);
    186183}
    187184
    188 
    189185bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
    190186{
    191 
    192     float dBin;
    193     int   nKeep, nSkip;
    194 
    195187    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    196188    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    … …  
    202194    //   we use an outlier rejection to avoid this bias
    203195
    204     // rflux = ten(0.4*fitMag);
    205     psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     196    // r2rflux = radius^2 * ten(0.4*fitMag);
     197    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    206198    for (int i = 0; i < psfTry->sources->n; i++) {
    207199        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
    208200            continue;
    209         rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
    210     }
    211 
    212     // find min and max of (1/flux):
    213     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    214     psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
    215 
    216     // build binned versions of rflux, metric
    217     dBin = (stats->max - stats->min) / 10.0;
    218     psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    219     psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    220     psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    221     psFree (stats);
    222 
    223     psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    224 
    225     // group data in daBin bins, measure lower 50% mean
    226     for (int i = 0; i < daBin->n; i++) {
    227 
    228         psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    229         tmp->n = 0;
    230 
    231         // accumulate data within bin range
    232         for (int j = 0; j < psfTry->sources->n; j++) {
    233             // masked for: bad model fit, outlier in parameters
    234             if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
    235                 continue;
    236 
    237             // skip points with extreme dA values
    238             if (fabs(psfTry->metric->data.F64[j]) > 0.5)
    239                 continue;
    240 
    241             // skip points outside of this bin
    242             if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    243                 continue;
    244             if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    245                 continue;
    246 
    247             tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
    248             tmp->n ++;
    249         }
    250 
    251         // is this a valid point?
    252         maskB->data.U8[i] = 0;
    253         if (tmp->n < 2) {
    254             maskB->data.U8[i] = 1;
    255             psFree (tmp);
    256             continue;
    257         }
    258 
    259         // dA values are contaminated with low outliers
    260         // measure statistics only on upper 50% of points
    261         // this would be easier if we could sort in reverse:
    262         //
    263         // psVectorSort (tmp, tmp);
    264         // tmp->n = 0.5*tmp->n;
    265         // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    266         // psVectorStats (stats, tmp, NULL, NULL, 0);
    267         // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    268 
    269         psVectorSort (tmp, tmp);
    270         nKeep = 0.5*tmp->n;
    271         nSkip = tmp->n - nKeep;
    272 
    273         psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    274         for (int j = 0; j < tmp2->n; j++) {
    275             tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    276         }
    277 
    278         stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    279         psVectorStats (stats, tmp2, NULL, NULL, 0);
    280         psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    281 
    282         daBin->data.F64[i] = stats->sampleMedian;
    283 
    284         psFree (stats);
    285         psFree (tmp);
    286         psFree (tmp2);
    287     }
    288 
    289     // linear fit to rfBin, daBin
    290     psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    291     psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    292     poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
    293 
    294     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    295     // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
    296 
    297     // XXX EAM : replace this when the above version is implemented
    298     // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
    299 
    300     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    301     psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    302 
    303     stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    304     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    305 
    306     psfTry->psf->ApResid = poly->coeff[0];
    307     psfTry->psf->dApResid = stats->clippedStdev;
    308     psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    309 
    310     psFree (rflux);
    311     psFree (rfBin);
    312     psFree (daBin);
    313     psFree (maskB);
    314     psFree (daBinFit);
    315     psFree (daResid);
    316     psFree (poly);
    317     psFree (stats);
    318     psFree (fitstat);
    319 
    320     return true;
    321 }
    322 
    323 bool pmPSFtryMetric_Alt (pmPSFtry *try
    324                          , float RADIUS)
    325 {
    326 
    327     // the measured (aperture - fit) magnitudes (dA == try->metric)
    328     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    329     //     dA = dAo + dsky/flux
    330     //   where flux is the flux of the star
    331     // we fit this trend to find the infinite flux aperture correction (dAo),
    332     //   the nominal sky bias (dsky), and the error on dAo
    333     // the values of dA are contaminated by stars with close neighbors in the aperture
    334     //   we use an outlier rejection to avoid this bias
    335 
    336     // rflux = ten(0.4*fitMag);
    337     psVector *rflux = psVectorAlloc (try
    338                                      ->sources->n, PS_TYPE_F64);
    339     for (int i = 0; i < try
    340                 ->sources->n; i++) {
    341             if (try
    342                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    343             rflux->data.F64[i] = pow(10.0, 0.4*try
    344                                      ->fitMag->data.F64[i]);
    345         }
    346 
    347     // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     201        r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
     202    }
     203
     204    // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
    348205    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    349 
    350     // XXX EAM
    351206    stats->min = 1.0;
    352207    stats->max = 3.0;
    353208    stats->clipIter = 3;
    354209
    355     // linear clipped fit to rfBin, daBin
    356     psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    357     poly = psVectorClipFitPolynomial1D (poly, stats, try
    358                                             ->mask, PSFTRY_MASK_ALL, try
    359                                                 ->metric, NULL, rflux);
    360     fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    361 
    362     try
    363         ->psf->ApResid = poly->coeff[0];
    364     try
    365         ->psf->dApResid = stats->sampleStdev;
    366     try
    367         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    368 
    369     fprintf (stderr, "*******************************************************************************\n");
    370 
    371     FILE *f;
    372     f = fopen ("apresid.dat", "w");
    373     if (f == NULL)
    374         psAbort ("pmPSFtry", "can't open output file");
    375 
    376     for (int i = 0; i < try
    377                 ->sources->n; i++) {
    378             fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
    379                          ->fitMag->data.F64[i], rflux->data.F64[i], try
    380                              ->metric->data.F64[i], try
    381                                  ->mask->data.U8[i]);
    382         }
    383     fclose (f);
    384 
    385     psFree (rflux);
     210    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
     211    // linear clipped fit of ApResid to r2rflux
     212    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     213    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux);
     214    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     215
     216    pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS);
     217    psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
     218    psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;
     219
     220    psfTry->psf->skyBias  = poly->coeff[1];
     221    psfTry->psf->ApResid  = poly->coeff[0];
     222    psfTry->psf->dApResid = stats->sampleStdev;
     223
     224    psFree (r2rflux);
    386225    psFree (poly);
    387226    psFree (stats);
    388227
    389     // psFree (daFit);
    390     // psFree (daResid);
    391 
    392228    return true;
    393229}
     230
     231/*
     232  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
     233  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
     234  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
     235*/
     236
Note: See TracChangeset for help on using the changeset viewer.