IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 24, 2007, 11:27:58 AM (19 years ago)
Author:
eugene
Message:

update from eam_branch_20070921

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14937 r15000  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-09-21 00:06:57 $
     7 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-09-24 21:27:58 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4444static void pmPSFtryFree (pmPSFtry *test)
    4545{
    46 
    47     if (test == NULL)
    48         return;
     46    if (test == NULL) return;
    4947
    5048    psFree (test->psf);
    … …  
    5856
    5957// allocate a pmPSFtry based on the desired sources and the model (identified by name)
    60 pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask)
     58pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options)
    6159{
    62 
    63     // validate the requested model name
    64     pmModelType type = pmModelClassGetType (modelName);
    65     if (type == -1) {
    66         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
    67         return NULL;
    68     }
    69 
    7060    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    71 
    72     test->psf       = pmPSFAlloc (type, poissonErrors, psfTrendMask);
     61    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
     62
     63    test->psf       = pmPSFAlloc (options);
    7364    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
    7465    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    7667    test->mask      = psVectorAlloc (sources->n, PS_TYPE_U8);
    7768
     69    psVectorInit (test->mask,        0);
     70    psVectorInit (test->metric,    0.0);
     71    psVectorInit (test->metricErr, 0.0);
     72    psVectorInit (test->fitMag,    0.0);
     73
    7874    test->sources   = psArrayAlloc (sources->n);
     75
    7976    for (int i = 0; i < sources->n; i++) {
    8077        test->sources->data[i] = pmSourceCopy (sources->data[i]);
    81         test->mask->data.U8[i]  = 0;
    82         test->metric->data.F32[i] = 0;
    83         test->metricErr->data.F32[i] = 0;
    84         test->fitMag->data.F32[i] = 0;
    85     }
    86 
    87     psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
     78    }
     79
    8880    return (test);
    8981}
    … …  
    9991
    10092// generate a pmPSFtry with a copy of the test PSF sources
    101 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark)
     93pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark)
    10294{
    10395    bool status;
    … …  
    10597    int Npsf = 0;
    10698
    107     pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask);
     99    // validate the requested model name
     100    options->type = pmModelClassGetType (modelName);
     101    if (options->type == -1) {
     102        psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
     103        return NULL;
     104    }
     105
     106    pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
    108107    if (psfTry == NULL) {
    109108        psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    … …  
    111110    }
    112111
    113     // stage 1:  fit an independent model (freeModel) to all sources
     112    // stage 1:  fit an EXT model to all candidates PSF sources
    114113    psTimerStart ("fit");
    115114    for (int i = 0; i < psfTry->sources->n; i++) {
    … …  
    123122        }
    124123
    125         // set object mask to define valid pixels
    126         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);
     124        // set object mask to define valid pixels -- XXX not unset?
     125        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
    127126
    128127        // fit model as EXT, not PSF
    … …  
    140139
    141140    // stage 2: construct a psf (pmPSF) from this collection of model fits
    142     if (!pmPSFFromPSFtry (psfTry, applyWeights)) {
     141    if (!pmPSFFromPSFtry (psfTry)) {
    143142        psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    144143        psFree(psfTry);
    … …  
    162161            continue;
    163162        }
    164         source->modelPSF->radiusFit = RADIUS;
    165 
    166         // set object mask to define valid pixels
    167         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);
     163        source->modelPSF->radiusFit = options->radius;
     164
     165        // set object mask to define valid pixels -- XXX not unset?
     166        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
    168167
    169168        // fit the PSF model to the source
    … …  
    241240
    242241    // XXX this function wants aperture radius for pmSourcePhotometry
    243     pmPSFtryMetric (psfTry, RADIUS);
     242    pmPSFtryMetric (psfTry, options->radius);
    244243    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
    245244              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    … …  
    339338/*****************************************************************************
    340339pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
    341 source->modelEXT entires.  The PSF ignores the first 4 (independent) model
     340source->modelEXT entries.  The PSF ignores the first 4 (independent) model
    342341parameters and constructs a polynomial fit to the remaining as a function of
    343342image coordinate.
    … …  
    345344Note: some of the array entries may be NULL (failed fits); ignore them.
    346345 *****************************************************************************/
    347 bool pmPSFFromPSFtry (pmPSFtry *psfTry, bool applyWeight)
     346bool pmPSFFromPSFtry (pmPSFtry *psfTry)
    348347{
    349348    pmPSF *psf = psfTry->psf;
    … …  
    355354
    356355    psVector *dz = NULL;
    357     if (applyWeight) {
     356    if (psf->poissonErrorsParams) {
    358357        dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    359358    }
    … …  
    377376    // more deviant than three sigma.
    378377    // mask is currently updated for each pass. this is inconsistent?
    379 
    380     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    381378
    382379    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    … …  
    389386    // threshold.
    390387
     388    // XXX this function needs to check the fit residuals and modify Nx,Ny as needed
     389
    391390    // convert the measured source shape paramters to polarization terms
    392391    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    … …  
    406405    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    407406    // This way, the parameters masked by one of the fits will be applied to the others
    408     for (int i = 0; i < stats->clipIter; i++) {
    409         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
    410         psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n);
    411         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
    412         psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n);
    413         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
    414         psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n);
    415     }
    416 
    417     // XXX temporary dump of the psf parameters
     407    for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     408
     409        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     410        pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz);
     411        psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
     412        pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz);
     413        psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
     414        pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz);
     415        psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     416    }
     417
     418    // test dump of the psf parameters
    418419    if (psTraceGetLevel("psModules.objects") >= 4) {
    419420        FILE *f = fopen ("pol.dat", "w");
     421        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    420422        for (int i = 0; i < e0->n; i++) {
    421             fprintf (f, "%f %f  :  %f %f %f  : %d\n",
     423            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    422424                     x->data.F32[i], y->data.F32[i],
    423                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]);
     425                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
     426                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
     427                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
     428                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
     429                     psfTry->mask->data.U8[i]);
    424430        }
    425431        fclose (f);
    … …  
    455461        // fit the collection of measured parameters to the PSF 2D model
    456462        // the mask is carried from previous steps and updated with this operation
    457         // the weight is either the flux error or NULL, depending on 'applyWeights'
    458         if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
     463        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
     464        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) {
    459465            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    460             psFree(stats);
    461466            psFree(x);
    462467            psFree(y);
    … …  
    485490
    486491            for (int i = 0; i < psf->params->n; i++) {
    487                 if (psf->params->data[i] == NULL)
    488                     continue;
     492                if (psf->params->data[i] == NULL) continue;
    489493                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
    490494            }
    … …  
    496500    }
    497501
    498     psFree (stats);
    499502    psFree (x);
    500503    psFree (y);
    501504    psFree (z);
    502505    psFree (dz);
    503     return (true);
     506    return true;
    504507}
    505508
Note: See TracChangeset for help on using the changeset viewer.