IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14980


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

updating to use pmTrend for the psf parameters

Location:
branches/eam_branch_20070921
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r14972 r14980  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.28.2.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-09-21 21:45:00 $
     8 *  @version $Revision: 1.28.2.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-09-22 02:20:11 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4545
    4646static void pmPSFOptionsFree (pmPSFOptions *options) {
     47
     48    if (!options) return;
     49
     50    psFree (options->image);
     51    psFree (options->stats);
    4752    return;
    4853}
    … …  
    5459
    5560    options->type          = 0;
    56     options->poissonErrors = false;
     61
     62    options->image         = NULL;
     63    options->stats         = NULL;
     64
    5765    options->psfTrendMode  = PM_TREND_NONE;
    5866    options->psfTrendNx    = 0;
    5967    options->psfTrendNy    = 0;
     68
     69    options->poissonErrors = false;
    6070}
    6171
    … …  
    136146    }
    137147    psf->params = psArrayAlloc(Nparams);
     148   
     149    // save the trend stats on the psf for use in pmPSFFromPSFtry
     150    psf->psfTrendStats = psMemIncrRefCounter (options->stats);
    138151
    139152    // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
    … …  
    144157    // PM_PAR_XPOS, etc)
    145158
    146     // XXX make this a user option
    147     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    148 
    149     // XXX define new structure to carry in the psf options
     159    // define the parameter trends
    150160    if (options->psfTrendMode != PM_TREND_NONE) {
    151161        for (int i = 0; i < psf->params->n; i++) {
    … …  
    155165            if (i == PM_PAR_YPOS) continue;
    156166
    157             // XXX need to set stats & image or else nXfine,nYfine
    158             pmTrend2D *param = pmTrend2DAlloc (options->psfTrendMode, image, options->psfTrendNx, options->psfTrendNx, stats);
    159             psf->params->data[i] = param;
     167            psf->params->data[i] = pmTrend2DAlloc (options->psfTrendMode, options->image, options->psfTrendNx, options->psfTrendNx, options->stats);
    160168        }
    161169    }
    162    
    163     psFree (stats);
    164170    return psf;
    165171}
    … …  
    287293
    288294// generate a psf model of the requested type, with fixed shape
     295// XXX update this.
    289296pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...)
    290297{

  • branches/eam_branch_20070921/psModules/src/objects/pmPSF.h

    r14969 r14980  
    66 * @author EAM, IfA
    77 *
    8  * @version $Revision: 1.16.2.1 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2007-09-21 18:55:12 $
     8 * @version $Revision: 1.16.2.2 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2007-09-22 02:20:11 $
    1010 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 */
    … …  
    3636    pmModelType type;                   ///< PSF Model in use
    3737    psArray *params;                    ///< Model parameters (psPolynomial2D)
     38    psStats *psfTrendStats;             ///< psf parameter trend clipping stats
    3839    psPolynomial1D *ChiTrend;           ///< Chisq vs flux fit (correction for systematic errors)
    3940    pmTrend2D *ApTrend;                 ///< ApResid vs (x,y)
    … …  
    5455typedef struct {
    5556    pmModelType   type;
    56     bool          poissonErrors;
     57    psImage      *image;                // image for which the PSF is defined
     58    psStats      *stats;                // psfTrend clipping stats
    5759    pmTrend2DMode psfTrendMode;
    5860    int           psfTrendNx;
    5961    int           psfTrendNy;
     62    bool          poissonErrors;
    6063} pmPSFOptions;
    6164
    … …  
    7578);
    7679
    77 bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option);
    78 
    79 pmPSFApTrendOptions pmPSFApTrendOptionFromName (char *name);
    80 
    8180double pmPSF_SXYfromModel (psF32 *modelPar);
    8281double pmPSF_SXYtoModel (psF32 *fittedPar);
    … …  
    8685
    8786bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes);
     87bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
     88
     89psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
    8890psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR);
    89 bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
    90 psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
    9191
    9292/// @}

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

    r14969 r14980  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.22.2.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-09-21 18:55:12 $
     8 *  @version $Revision: 1.22.2.2 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-09-22 02:20:11 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    148148
    149149// for a pmPSF supplied on the analysis metadata, we write out
     150// if needed:
     151//   - a PHU blank header
    150152// - image header        : FITS Image NAXIS = 0
    151 // - psf table (+header) : FITS Table
    152 // XXX if we are using psImageMap, this is also a FITS Image, with header data
    153 // - psf resid (+header) : FITS Image
    154 // if needed, we also write out a PHU blank header
     153// if (trendMode == MAP)
     154//   - psf resid (+header) : FITS Image
     155// else
     156//   - psf table (+header) : FITS Table
    155157bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    156158{
    … …  
    216218
    217219    // write out the IMAGE header segment
    218     // this header block is new, write it to disk
     220    // if this header block is new, write it to disk
    219221    if (hdu->header != file->header) {
    220222        // add EXTNAME, EXTHEAD, EXTTYPE to header
    … …  
    244246
    245247    // write the PSF model parameters in a FITS table
     248    assert (psf->psfTrendMode != PM_TREND_NONE);
    246249    {
    247250        // we need to write a header for the table,
    … …  
    267270        }
    268271
     272        char *modeName = pmTrend2DModeToString(trend->mode);
     273
    269274        // other required information describing the PSF
    270         psMetadataAddF32  (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
    271         psMetadataAddF32  (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
    272         psMetadataAddF32  (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
    273         psMetadataAddS32  (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
     275        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
     276        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
     277        psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
     278        psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
     279        psMetadataAddStr (header, PS_LIST_TAIL, "MODE",     0, "", modeName);
     280        psFree (modeName);
    274281
    275282        // XXX can we drop this now?
    … …  
    279286        psArray *psfTable = psArrayAllocEmpty (100);
    280287        for (int i = 0; i < nPar; i++) {
    281             psPolynomial2D *poly = psf->params->data[i];
    282             if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
    283             for (int ix = 0; ix <= poly->nX; ix++) {
    284                 for (int iy = 0; iy <= poly->nY; iy++) {
    285 
    286                     psMetadata *row = psMetadataAlloc ();
    287                     psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
    288                     psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
    289                     psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
    290                     psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
    291                     psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
    292                     psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
    293 
    294                     psArrayAdd (psfTable, 100, row);
    295                     psFree (row);
     288            pmTrend2D *trend = psf->params->data[i];
     289            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
     290
     291            if (trend->mode == PM_TREND_MAP) {
     292                // write the image components into a table: this is needed because they may each be a different size
     293                psImageMap *map = trend->map;
     294                for (int ix = 0; ix <= map->map->numCols; ix++) {
     295                    for (int iy = 0; iy <= map->map->numRows; iy++) {
     296                        psMetadata *row = psMetadataAlloc ();
     297                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     298                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     299                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     300                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
     301                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[ix][iy]);
     302                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
     303
     304                        psArrayAdd (psfTable, 100, row);
     305                        psFree (row);
     306                    }
     307                }
     308            } else {
     309                // write the polynomial components into a table
     310                psPolynomial2D *poly = trend->poly;
     311                for (int ix = 0; ix <= poly->nX; ix++) {
     312                    for (int iy = 0; iy <= poly->nY; iy++) {
     313                        psMetadata *row = psMetadataAlloc ();
     314                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     315                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     316                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     317                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
     318                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
     319                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
     320
     321                        psArrayAdd (psfTable, 100, row);
     322                        psFree (row);
     323                    }
    296324                }
    297325            }
    … …  
    535563        int nYorder = psMetadataLookupS32 (&status, header, name);
    536564        if (!status) {
    537             psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i);
     565            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
    538566            return false;
    539567        }
     568        snprintf (name, 9, "PAR%02d_MD", i);
     569        char *modeName = psMetadataLookupStr (&status, header, name);
     570        if (!status) {
     571            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
     572            return false;
     573        }
     574        // XXX allocate pmTrend2d based on the values here
    540575        psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
    541576    }
    … …  
    549584    // XXX can we drop this now?
    550585    psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");
     586
     587    char *modeName = psMetadataLookupStr (&status, header, "MODE");
     588    trend->mode = pmTrend2DModeToString(modeName);
    551589
    552590    // read the raw table data
    … …  
    561599        // XXX sanity check here
    562600
    563         psPolynomial2D *poly = psf->params->data[iPar];
    564         if (poly == NULL) {
    565             psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar);
     601        pmTrend2D *trend = psf->params->data[iPar];
     602        if (trend == NULL) {
     603            psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
    566604            return false;
    567605        }
    568606
    569         poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    570         poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
    571         poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     607        if (trend->mode == PM_TREND_MAP) {
     608            psImageMap *map = trend->map;
     609            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     610            map->error->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "ERROR");
     611            // poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     612        } else {
     613            psPolynomial2D *poly = trend->poly;
     614            poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     615            poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
     616            poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     617        }
    572618    }
    573619    psFree (header);

  • 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            }

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

    r14969 r14980  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2007-09-21 18:55:12 $
     5 *  @version $Revision: 1.2.2.2 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2007-09-22 02:20:11 $
    77 *
    88 *  Copyright 2004 Institute for Astronomy, University of Hawaii
    … …  
    179179    return result;
    180180}
     181
     182psString pmTrend2DModeToString (pmTrend2DMode mode) {
     183
     184    switch (mode) {
     185      case PM_TREND_NONE:
     186        psString name = psStringCopy ("NONE");
     187        break;
     188      case PM_TREND_POLY_ORD:
     189        psString name = psStringCopy ("POLY_ORD");
     190        break;
     191      case PM_TREND_POLY_CHEB:
     192        psString name = psStringCopy ("POLY_CHEB");
     193        break;
     194      case PM_TREND_MAP:
     195        psString name = psStringCopy ("MAP");
     196        break;
     197      default:
     198        psAbort ("invalid mode %d", mode);
     199    }
     200    return name;
     201}
     202
     203pmTrend2DMode pmTrend2DModeFromString (psString name) {
     204
     205    if (!name) return PM_TREND_NONE;
     206
     207    if (!strcasecmp (name, "NONE")) {
     208        return PM_TREND_NONE;
     209    }
     210    if (!strcasecmp (name, "POLY_ORD")) {
     211        return PM_TREND_POLY_ORD;
     212    }
     213    if (!strcasecmp (name, "POLY_CHEB")) {
     214        return PM_TREND_POLY_CHEB;
     215    }
     216    if (!strcasecmp (name, "MAP")) {
     217        return PM_TREND_MAP;
     218    }
     219    psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name);
     220    return PM_TREND_NONE;
     221}

  • branches/eam_branch_20070921/psModules/src/objects/pmTrend2D.h

    r14972 r14980  
    55 * @author EAM, IfA
    66 *
    7  * @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2007-09-21 21:45:00 $
     7 * @version $Revision: 1.2.2.2 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2007-09-22 02:20:11 $
    99 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1010 */
    … …  
    4242psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y);
    4343
     44psString pmTrend2DModeToString (pmTrend2DMode mode);
     45pmTrend2DMode pmTrend2DModeFromString (psString name);
     46
    4447/// @}
    4548# endif

  • branches/eam_branch_20070921/psphot/src/psphotChoosePSF.c

    r14971 r14980  
    3535
    3636    // how to model the PSF variations across the field
    37     char *psfTrendMode = psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE");
    38     assert (status);
    39     options->psfTrendMode = pmTrend2DModeFromString (psfTrendMode);
     37    options->psfTrendMode = pmTrend2DModeFromString (psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE"));
     38    assert (status);
     39   
     40    options->psfTrendNx = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NX");
     41    assert (status);
     42
     43    options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY");
     44    assert (status);
    4045
    4146    // get the fixed PSF fit radius
    … …  
    4449    assert (status);
    4550
    46     pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), POISSON_ERRORS);
    47 
    48     psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
    49     if (!psfTrendMask) {
    50         psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
    51         return NULL;
    52     }
     51    options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     52    options->image = psMemIncrRefCounter (readout->image);
     53
     54    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotom);
    5355
    5456    psArray *stars = psArrayAllocEmpty (sources->n);
    … …  
    7173        psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
    7274        psFree (stars);
    73         psFree (psfTrendMask);
    7475        return NULL;
    7576    }
    … …  
    9697        psMetadataItem *item = psListGetAndIncrement (iter);
    9798        char *modelName = item->data.V;
    98         // models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
    9999        models->data[i] = pmPSFtryModel (stars, modelName, options, maskVal, mark);
    100100    }
    … …  
    103103    psFree (list);
    104104    psFree (stars);
    105     psFree (psfTrendMask);
    106105
    107106    // select the best of the models
    … …  
    134133    if (psTraceGetLevel("psphot") >= 5) {
    135134        for (int i = PM_PAR_SXX; i < try->psf->params->n; i++) {
    136             psPolynomial2D *poly = try->psf->params->data[i];
    137             for (int nx = 0; nx <= poly->nX; nx++) {
    138                 for (int ny = 0; ny <= poly->nY; ny++) {
    139                     if (poly->mask[nx][ny]) continue;
    140                     fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
    141                 }
    142             }
     135            // XXX re-write the output or some other status info
    143136        }
    144137    }
    … …  
    289282        pmSourcesWritePSFs (try->sources, "psfstars.dat");
    290283        pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
    291         psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
    292         psMetadataConfigWrite (psfData, "psfmodel.dat");
    293         psFree (psfData);
     284        // XXX need alternative output function
     285        // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
     286        // psMetadataConfigWrite (psfData, "psfmodel.dat");
    294287        psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
    295288        exit (0);
Note: See TracChangeset for help on using the changeset viewer.