IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 21, 2007, 4:20:17 PM (19 years ago)
Author:
eugene
Message:

updating to use pmTrend for the psf parameters

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14972 r14980  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.45.2.1 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-09-21 21:45:00 $
     7 *  @version $Revision: 1.45.2.2 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-09-22 02:20:11 $
    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, pmPSFOptions *options)
     58pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options)
    6159{
    62 
    63     // validate the requested model name
    64     options->type = pmModelClassGetType (modelName);
    65     if (options->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));
     61    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    7162
    7263    test->psf       = pmPSFAlloc (options);
    … …  
    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}
    … …  
    10597    int Npsf = 0;
    10698
    107     // XXX use the options structure here instead?
    108 
    109     pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, options);
     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);
    110107    if (psfTry == NULL) {
    111108        psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    … …  
    113110    }
    114111
    115     // stage 1:  fit an independent model (freeModel) to all sources
     112    // stage 1:  fit an EXT model to all candidates PSF sources
    116113    psTimerStart ("fit");
    117114    for (int i = 0; i < psfTry->sources->n; i++) {
    … …  
    125122        }
    126123
    127         // set object mask to define valid pixels
     124        // set object mask to define valid pixels -- XXX not unset?
    128125        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
    129126
    … …  
    142139
    143140    // stage 2: construct a psf (pmPSF) from this collection of model fits
     141    // XXX take 'applyWeights' from the psf options?
    144142    if (!pmPSFFromPSFtry (psfTry, options->applyWeights)) {
    145143        psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    … …  
    166164        source->modelPSF->radiusFit = options->radius;
    167165
    168         // set object mask to define valid pixels
     166        // set object mask to define valid pixels -- XXX not unset?
    169167        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
    170168
    … …  
    380378    // mask is currently updated for each pass. this is inconsistent?
    381379
    382     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    383 
    384380    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    385381    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
    … …  
    391387    // threshold.
    392388
     389    // XXX this function needs to check the fit residuals and modify Nx,Ny as needed
     390
    393391    // convert the measured source shape paramters to polarization terms
    394392    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    … …  
    408406    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    409407    // This way, the parameters masked by one of the fits will be applied to the others
    410     for (int i = 0; i < stats->clipIter; i++) {
    411         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
    412         psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n);
    413         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
    414         psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n);
    415         psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
    416         psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n);
    417     }
    418 
    419     // XXX temporary dump of the psf parameters
     408    for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     409
     410        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     411        pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz);
     412        psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
     413        pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz);
     414        psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
     415        pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz);
     416        psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
     417    }
     418
     419    // test dump of the psf parameters
    420420    if (psTraceGetLevel("psModules.objects") >= 4) {
    421421        FILE *f = fopen ("pol.dat", "w");
     422        fprintf (f, "x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n",
    422423        for (int i = 0; i < e0->n; i++) {
    423             fprintf (f, "%f %f  :  %f %f %f  : %d\n",
     424            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    424425                     x->data.F32[i], y->data.F32[i],
    425                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]);
     426                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
     427                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
     428                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
     429                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
     430                     psfTry->mask->data.U8[i]);
    426431        }
    427432        fclose (f);
    … …  
    458463        // the mask is carried from previous steps and updated with this operation
    459464        // the weight is either the flux error or NULL, depending on 'applyWeights'
    460         if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
     465        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) {
    461466            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    462467            psFree(stats);
    … …  
    487492
    488493            for (int i = 0; i < psf->params->n; i++) {
    489                 if (psf->params->data[i] == NULL)
    490                     continue;
     494                if (psf->params->data[i] == NULL) continue;
    491495                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
    492496            }
Note: See TracChangeset for help on using the changeset viewer.