IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9886


Ignore:
Timestamp:
Nov 6, 2006, 11:14:00 PM (20 years ago)
Author:
magnier
Message:

convert apResid measurement to use pmSourceMagnitudes API for consistency; now using a new source for each pmPSFtry; dropped modelEXT and modelPSF arrays; instead use the source->modelEXT,PSF elements of the pmPSFtry

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9841 r9886  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-11-03 15:12:24 $
     7 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-11-07 09:14:00 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4646    psFree (test->psf);
    4747    psFree (test->sources);
    48     psFree (test->modelEXT);
    49     psFree (test->modelPSF);
     48    // psFree (test->modelEXT);
     49    // psFree (test->modelPSF);
    5050    psFree (test->metric);
    5151    psFree (test->metricErr);
    … …  
    6969
    7070    test->psf       = pmPSFAlloc (type, poissonErrors, psfTrendMask);
    71     test->sources   = psMemIncrRefCounter(sources);
    72     test->modelEXT  = psArrayAlloc (sources->n);
    73     test->modelPSF  = psArrayAlloc (sources->n);
    7471    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F64);
    7572    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F64);
    7673    test->fitMag    = psVectorAlloc (sources->n, PS_TYPE_F64);
    7774    test->mask      = psVectorAlloc (sources->n, PS_TYPE_U8);
    78 
    79     for (int i = 0; i < test->modelEXT->n; i++) {
     75    // test->modelEXT  = psArrayAlloc (sources->n);
     76    // test->modelPSF  = psArrayAlloc (sources->n);
     77
     78    test->sources   = psArrayAlloc (sources->n);
     79    for (int i = 0; i < sources->n; i++) {
     80        test->sources->data[i] = pmSourceCopy (sources->data[i]);
    8081        test->mask->data.U8[i]  = 0;
    81         test->modelEXT->data[i] = NULL;
    82         test->modelPSF->data[i] = NULL;
    8382        test->metric->data.F64[i] = 0;
    8483        test->metricErr->data.F64[i] = 0;
    8584        test->fitMag->data.F64[i] = 0;
     85        // test->modelEXT->data[i] = NULL;
     86        // test->modelPSF->data[i] = NULL;
    8687    }
    8788
    … …  
    99100// mask values indicate the reason the source was rejected:
    100101
     102// generate a pmPSFtry with a copy of the test PSF sources
    101103pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask)
    102104{
    103105    bool status;
    104     float obsMag;
    105     float fitMag;
    106106    float x;
    107107    float y;
    … …  
    111111    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask);
    112112    if (psfTry == NULL) {
    113         psError (PS_ERR_UNKNOWN, false, "failed to create allocate psf model");
     113        psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    114114        return NULL;
    115115    }
    … …  
    120120
    121121        pmSource *source = psfTry->sources->data[i];
    122         pmModel  *model = pmSourceModelGuess (source, psfTry->psf->type);
    123         if (model == NULL) {
     122        source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
     123        if (source->modelEXT == NULL) {
    124124            psError(PS_ERR_UNKNOWN, false, "failed to build model");
    125125            return NULL;
    … …  
    132132
    133133        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
    134         status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT);
     134        status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT);
    135135        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    136136
    … …  
    138138        if (!status) {
    139139            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
    140             psFree (model);
    141             continue;
    142         }
    143         psfTry->modelEXT->data[i] = model;
     140            continue;
     141        }
     142        // XXX is this still needed?
     143        // psfTry->modelEXT->data[i] = psMemIncrRefCounter (source->modelEXT);
    144144        Next ++;
    145145    }
    … …  
    148148
    149149    // stage 2: construct a psf (pmPSF) from this collection of model fits
    150     pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);
     150    pmPSFFromPSFtry (psfTry);
    151151
    152152    // stage 3: refit with fixed shape parameters
    … …  
    158158
    159159        pmSource *source = psfTry->sources->data[i];
    160         pmModel  *modelEXT = psfTry->modelEXT->data[i];
    161160
    162161        // set shape for this model based on PSF
    163         pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
     162        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
     163        source->modelPSF->radiusFit = RADIUS;
    164164        x = source->peak->x;
    165165        y = source->peak->y;
    166166
     167        // set the mask and fit the PSF model
    167168        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
    168         status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF);
     169        status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF);
     170        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    169171
    170172        // skip poor fits
    171173        if (!status) {
    172174            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
    173             psFree (modelPSF);
    174             goto next_source;
     175            continue;
    175176        }
    176177
    177178        // otherwise, save the resulting model
    178         psfTry->modelPSF->data[i] = modelPSF;
    179 
    180         // I'm not using the pmSourceMagnitudes API, why?
     179        // psfTry->modelPSF->data[i] = psMemIncrRefCounter (source->modelPSF);
     180
    181181        // XXX : use a different aperture radius from the fit radius?
    182         // XXX : use a different estimator for the local sky?
    183         // XXX : pass 'source' as input?
    184         if (!pmSourcePhotometryModel (&fitMag, modelPSF)) {
     182        if (!pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP)) {
    185183            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    186             goto next_source;
    187         }
    188         if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) {
    189             psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    190             goto next_source;
    191         }
    192 
    193         psfTry->metric->data.F64[i] = obsMag - fitMag;
    194         psfTry->metricErr->data.F64[i] = modelPSF->dparams->data.F32[PM_PAR_I0] / modelPSF->params->data.F32[PM_PAR_I0];
    195 
    196         psfTry->fitMag->data.F64[i] = fitMag;
     184            continue;
     185        }
     186
     187        psfTry->fitMag->data.F64[i] = source->psfMag;
     188        psfTry->metric->data.F64[i] = source->apMag - source->psfMag;
     189        psfTry->metricErr->data.F64[i] = source->errMag;
    197190        Npsf ++;
    198 
    199 next_source:
    200         psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    201 
    202191    }
    203192    psfTry->psf->nPSFstars = Npsf;
    … …  
    212201    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK);
    213202
    214     // write sources with models first
     203    // generate the x and y vectors, and mask missing models
    215204    for (int i = 0; i < psfTry->sources->n; i++) {
    216         pmModel *model = psfTry->modelPSF->data[i];
    217         if (model == NULL) {
     205        pmSource *source = psfTry->sources->data[i];
     206        if (source->modelPSF == NULL) {
    218207            flux->data.F64[i] = 0.0;
    219208            chisq->data.F64[i] = 0.0;
    220209            mask->data.U8[i] = 1;
    221210        } else {
    222             flux->data.F64[i] = model->params->data.F32[PM_PAR_I0];
    223             chisq->data.F64[i] = model[0].chisq/model[0].nDOF;
     211            flux->data.F64[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
     212            chisq->data.F64[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
    224213            mask->data.U8[i] = 0;
    225214        }
    226215    }
    227     // use 3hi/3lo sigma clipping on the r2rflux vs metric fit
     216
     217    // use 3hi/3lo sigma clipping on the chisq fit
    228218    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    229219    stats->clipSigma = 3.0;
    230220    stats->clipIter = 3;
    231     // linear clipped fit of ApResid to r2rflux
    232     psfTry->psf->ChiTrend =
    233         psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
     221
     222    // linear clipped fit of chisq trend vs flux
     223    psfTry->psf->ChiTrend = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
    234224    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    235225
    … …  
    246236
    247237    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
    248         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]*pow(10000,i));
     238        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
     239                  psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
     240                  psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
    249241    }
    250242
    … …  
    342334*/
    343335
     336/*****************************************************************************
     337pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
     338source->modelEXT entires.  The PSF ignores the first 4 (independent) model
     339parameters and constructs a polynomial fit to the remaining as a function of
     340image coordinate.
     341    input: psfTry with fitted source->modelEXT collection, pre-allocated psf
     342Note: some of the array entries may be NULL (failed fits); ignore them.
     343 *****************************************************************************/
     344bool pmPSFFromPSFtry (pmPSFtry *psfTry)
     345{
     346    pmPSF *psf = psfTry->psf;
     347
     348    // construct the fit vectors from the collection of objects
     349    psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     350    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     351    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     352    psVector *dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     353
     354    // construct the x,y terms
     355    for (int i = 0; i < psfTry->sources->n; i++) {
     356        pmSource *source = psfTry->sources->data[i];
     357        if (source->modelEXT == NULL)
     358            continue;
     359
     360        // use F64 for polynomial fitting
     361        x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
     362        y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
     363    }
     364
     365    // we are doing a robust fit.  after each pass, we drop points which are
     366    // more deviant than three sigma.
     367    // mask is currently updated for each pass. this is inconsistent?
     368
     369    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     370
     371    // skip the unfitted parameters (X, Y, Io, Sky)
     372    for (int i = 0; i < psf->params_NEW->n; i++) {
     373        if (i == PM_PAR_SKY)
     374            continue;
     375        if (i == PM_PAR_I0)
     376            continue;
     377        if (i == PM_PAR_XPOS)
     378            continue;
     379        if (i == PM_PAR_YPOS)
     380            continue;
     381
     382        // select the per-object fitted data for this parameter
     383        for (int j = 0; j < psfTry->sources->n; j++) {
     384            pmSource *source = psfTry->sources->data[i];
     385            if (source->modelEXT == NULL)
     386                continue;
     387
     388            z->data.F64[j] = source->modelEXT->params->data.F32[i];
     389            dz->data.F64[j] = 1; // use the model-fitted error? or S/N?
     390
     391            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
     392            if (i == PM_PAR_SXY) {
     393                z->data.F64[j] = pmPSF_SXYfromModel (source->modelEXT->params->data.F32);
     394            }
     395        }
     396        psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, dz, x, y);
     397        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, psfTry->mask->n, i);
     398    }
     399
     400    // XXX test dump of star parameters vs position (compare with fitted values)
     401    if (0) {
     402        FILE *f = fopen ("params.dat", "w");
     403
     404        for (int j = 0; j < psfTry->sources->n; j++) {
     405            pmSource *source = psfTry->sources->data[j];
     406            if (source == NULL)
     407                continue;
     408
     409            pmModel *model = source->modelEXT;
     410            if (model == NULL)
     411                continue;
     412
     413            pmModel *modelPSF = pmModelFromPSF (model, psf);
     414
     415            fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     416
     417            for (int i = 0; i < psf->params_NEW->n; i++) {
     418                if (psf->params_NEW->data[i] == NULL)
     419                    continue;
     420                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
     421            }
     422            fprintf (f, "%f %d\n", model->chisq, model->nIter);
     423        }
     424        fclose (f);
     425    }
     426
     427    psFree (stats);
     428    psFree (x);
     429    psFree (y);
     430    psFree (z);
     431    psFree (dz);
     432    return (true);
     433}
     434
Note: See TracChangeset for help on using the changeset viewer.