IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 11176


Ignore:
Timestamp:
Jan 18, 2007, 6:56:07 PM (19 years ago)
Author:
eugene
Message:

convert to PM_PAR_XXX names, added running count of each type of fit attempted

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceFits.c

    r10268 r11176  
    33// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
    44// XXX this function does not call pmSourceFitModelInit : fix this?
     5
     6static int NfitPSF = 0;
     7static int NfitBlend = 0;
     8static int NfitDBL = 0;
     9static int NfitEXT = 0;
     10
     11bool psphotFitInit () {
     12    psTimerStart ("psphot.fits");
     13    return true;
     14}
     15
     16bool psphotFitSummary () {
     17
     18    fprintf (stderr, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
     19             NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
     20    return true;
     21}
    522
    623bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) {
    … …  
    1431    }
    1532    psTrace ("psphot", 5, "trying blend...\n");
     33    return false;
    1634
    1735    // save the PSF model from the Ensemble fit
    1836    pmModel *PSF = pmModelCopy (source->modelPSF);
    19     if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary");
    20 
    21     x = PSF->params->data.F32[2];
    22     y = PSF->params->data.F32[3];
     37
     38    if (isnan(PSF->params->data.F32[PM_PAR_I0])) psAbort ("psphot", "nan in blend fit primary");
     39
     40    x = PSF->params->data.F32[PM_PAR_XPOS];
     41    y = PSF->params->data.F32[PM_PAR_YPOS];
    2342
    2443    psArray *modelSet = psArrayAllocEmpty (source->blends->n + 1);
    … …  
    3453
    3554        // find the blend which is furthest from source
    36         dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
     55        dR = PS_MAX (dR, hypot (blend->peak->xf - x, blend->peak->yf - y));
    3756
    3857        // create the model and guess parameters for this blend
    … …  
    4463
    4564        // XXX assume local sky is 0.0?
    46         model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
    47         if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
    48         model->params->data.F32[2] = blend->peak->x;
    49         model->params->data.F32[3] = blend->peak->y;
     65        model->params->data.F32[PM_PAR_I0] = blend->peak->flux;
     66        model->params->data.F32[PM_PAR_XPOS] = blend->peak->xf;
     67        model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf;
     68
     69        // these should never be invalid values
     70        // XXX drop these tests eventually
     71        if (isnan(model->params->data.F32[PM_PAR_I0])) psAbort ("psphot", "nan in blend fit");
     72        if (isnan(model->params->data.F32[PM_PAR_XPOS])) psAbort ("psphot", "nan in blend fit");
     73        if (isnan(model->params->data.F32[PM_PAR_YPOS])) psAbort ("psphot", "nan in blend fit");
    5074
    5175        // add this blend to the list
    … …  
    6690
    6791    // correct model chisq for flux trend
    68     double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     92    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    6993    PSF->chisqNorm = PSF->chisq / chiTrend;
    7094
    … …  
    7599
    76100        // correct model chisq for flux trend
    77         chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
     101        chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[PM_PAR_I0]);
    78102        model->chisqNorm = model->chisq / chiTrend;
    79103
    … …  
    92116        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
    93117    }
     118    NfitBlend += modelSet->n;
     119
    94120    // evaluate the primary object
    95121    if (!psphotEvalPSF (source, PSF)) {
    … …  
    117143    double chiTrend;
    118144
     145    NfitPSF ++;
     146
    119147    // save the PSF model from the Ensemble fit
    120148    pmModel *PSF = pmModelCopy (source->modelPSF);
    … …  
    124152    psphotCheckRadiusPSF (readout, source, PSF);
    125153
    126     x = PSF->params->data.F32[2];
    127     y = PSF->params->data.F32[3];
     154    x = PSF->params->data.F32[PM_PAR_XPOS];
     155    y = PSF->params->data.F32[PM_PAR_YPOS];
    128156
    129157    // fit PSF model (set/unset the pixel mask)
    … …  
    133161
    134162    // correct model chisq for flux trend
    135     chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     163    chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    136164    PSF->chisqNorm = PSF->chisq / chiTrend;
    137165
    … …  
    239267    psFree (DBL);
    240268    pmModelSub (source->pixels, source->mask, EXT, false, false);
    241     psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[2], EXT->params->data.F32[3]);
     269    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
    242270
    243271    // save new model
    … …  
    252280    pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false);
    253281    pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false);
    254     psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]);
     282    psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
    255283
    256284    // drop old model, save new second model...
    … …  
    287315    psArray *modelSet;
    288316
     317    NfitDBL ++;
     318
    289319    // make a guess at the position of the two sources
    290320    moments.x2 = source->moments->Sx;
    … …  
    305335
    306336    DBL = pmModelCopy (PSF);
    307     DBL->params->data.F32[1] *= 0.5;
    308     DBL->params->data.F32[2] = source->peak->xf + dx;
    309     DBL->params->data.F32[3] = source->peak->yf + dy;
     337    DBL->params->data.F32[PM_PAR_I0] *= 0.5;
     338    DBL->params->data.F32[PM_PAR_XPOS] = source->peak->xf + dx;
     339    DBL->params->data.F32[PM_PAR_YPOS] = source->peak->yf + dy;
    310340    modelSet->data[0] = DBL;
    311341
    312342    DBL = pmModelCopy (PSF);
    313     DBL->params->data.F32[1] *= 0.5;
    314     DBL->params->data.F32[2] = source->peak->xf - dx;
    315     DBL->params->data.F32[3] = source->peak->yf - dy;
     343    DBL->params->data.F32[PM_PAR_I0] *= 0.5;
     344    DBL->params->data.F32[PM_PAR_XPOS] = source->peak->xf - dx;
     345    DBL->params->data.F32[PM_PAR_YPOS] = source->peak->yf - dy;
    316346    modelSet->data[1] = DBL;
    317347
    … …  
    331361    float x, y;
    332362
     363    NfitEXT ++;
     364
    333365    // use the source moments, etc to guess basic model parameters
    334366    pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT);
    … …  
    337369    psphotCheckRadiusEXT (readout, source, EXT);
    338370
    339     x = EXT->params->data.F32[2];
    340     y = EXT->params->data.F32[3];
     371    x = EXT->params->data.F32[PM_PAR_XPOS];
     372    y = EXT->params->data.F32[PM_PAR_YPOS];
    341373
    342374    if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) {
Note: See TracChangeset for help on using the changeset viewer.