IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 7, 2005, 6:32:40 PM (21 years ago)
Author:
eugene
Message:

cleanup, organization

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/pspsf.c

    r4954 r4977  
    11# include "psphot.h"
    22
    3 static void pmPSFFree (pmPSF *psf) {
    4 
    5   if (psf == NULL) return;
    6 
    7   psFree (psf->params);
    8   return;
    9 }
    10 
    11 // a PSF always has 4 parameters fewer than the equivalent model
    12 pmPSF *pmPSFAlloc (pmModelType type) {
    13 
    14     int Nparams;
    15 
    16     pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
    17 
    18     psf->type = type;
    19     psf->chisq = 0.0;
    20 
    21     Nparams = pmModelParameterCount (type);
    22     if (!Nparams) {
    23         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    24         return(NULL);
    25     }     
    26 
    27     psf->params = psArrayAlloc (Nparams - 4);
    28     for (int i = 0; i < psf->params->n; i++) {
    29         // we need a way to set the allowed range of orders
    30         psf->params->data[i] = Polynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
    31     }
    32 
    33     psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
    34     return(psf);
    35 }
    36 
    37 static void pmPSF_TestFree (pmPSF_Test *test) {
     3// ********  pmPSFtry functions  **************************************************
     4// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
     5// * fitted version of the model, the pmPSF fit to the fitted model parameters,
     6// * and the PSF fits to the source. It also includes the statistics from the
     7// * fits, both the individual sources, and the collection
     8
     9// free a pmPSFtry structure
     10static void pmPSFtryFree (pmPSFtry *test) {
    3811
    3912  if (test == NULL) return;
    … …  
    4922}
    5023
    51 pmPSF_Test *pmPSF_TestAlloc (psArray *sources, char *modelName) {
    52 
    53     pmPSF_Test *test = (pmPSF_Test *) psAlloc(sizeof(pmPSF_Test));
     24// allocate a pmPSFtry based on the desired sources and the model (identified by name)
     25pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) {
     26
     27    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    5428
    5529    // XXX probably need to increment ref counter
    … …  
    7448    }   
    7549
    76     psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree);
     50    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    7751    return (test);
    7852}
    7953
    80 // test->mask values indicate reason source was rejected:
     54// build a pmPSFtry for the given model:
     55// - fit each source with the free-floating model
     56// - construct the pmPSF from the collection of models
     57// - fit each source with the PSF-parameter models
     58// - measure the pmPSF quality metric (dApResid)
     59
     60// sources used in for pmPSFtry may be masked by the analysis
     61// mask values indicate the reason the source was rejected:
    8162// 1: outlier in psf polynomial fit
    8263// 2: flt model failed to converge
    8364// 3: psf model failed to converge
    8465// 4: invalid source photometry
    85 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS)
     66// XXX EAM : use an enum for these values?
     67
     68pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
    8669{
    8770    bool status;
    … …  
    9376    int Npsf = 0;
    9477
    95     pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
     78    pmPSFtry *try = pmPSFtryAlloc (sources, modelName);
    9679
    9780    // stage 1:  fit an independent model (freeModel) to all sources
    9881    psTimerStart ("fit");
    99     for (int i = 0; i < test->sources->n; i++) {
    100 
    101         pmSource *source = test->sources->data[i];
    102         pmModel  *model  = pmSourceModelGuess (source, test->modelType);
     82    for (int i = 0; i < try->sources->n; i++) {
     83
     84        pmSource *source = try->sources->data[i];
     85        pmModel  *model  = pmSourceModelGuess (source, try->modelType);
    10386        x = source->peak->x;
    10487        y = source->peak->y;
    … …  
    11295        // exclude the poor fits
    11396        if (!status) {
    114           test->mask->data.U8[i] = 2;
     97          try->mask->data.U8[i] = 2;
    11598          psFree (model);
    11699          continue;
    117100        }
    118         test->modelFLT->data[i] = model;
     101        try->modelFLT->data[i] = model;
    119102        Nflt ++;
    120103    }
    121     psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    122     psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
     104    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     105    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    123106
    124107    // make this optional?
    125     // DumpModelFits (test->modelFLT, "modelsFLT.dat");
     108    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
    126109
    127110    // stage 2: construct a psf (pmPSF) from this collection of model fits
    128     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
     111    pmPSFFromModels (try->psf, try->modelFLT, try->mask);
    129112
    130113    // stage 3: refit with fixed shape parameters
    131114    psTimerStart ("fit");
    132     for (int i = 0; i < test->sources->n; i++) {
     115    for (int i = 0; i < try->sources->n; i++) {
    133116        // masked for: bad model fit, outlier in parameters
    134         if (test->mask->data.U8[i]) continue;
    135 
    136         pmSource *source = test->sources->data[i];
    137         pmModel  *modelFLT = test->modelFLT->data[i];
     117        if (try->mask->data.U8[i]) continue;
     118
     119        pmSource *source = try->sources->data[i];
     120        pmModel  *modelFLT = try->modelFLT->data[i];
    138121
    139122        // set shape for this model based on PSF
    140         pmModel *modelPSF = pmModelFromPSF (modelFLT, test->psf);
     123        pmModel *modelPSF = pmModelFromPSF (modelFLT, try->psf);
    141124        x = source->peak->x;
    142125        y = source->peak->y;
    … …  
    147130        // skip poor fits
    148131        if (!status) {
    149             test->mask->data.U8[i] = 3;
     132            try->mask->data.U8[i] = 3;
    150133            psFree (modelPSF);
    151134            goto next_source;
    … …  
    153136
    154137        // otherwise, save the resulting model
    155         test->modelPSF->data[i] = modelPSF;
     138        try->modelPSF->data[i] = modelPSF;
    156139
    157140        // XXX : use a different aperture radius from the fit radius?
    … …  
    159142        // XXX : pass 'source' as input?
    160143        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
    161             test->mask->data.U8[i] = 4;
     144            try->mask->data.U8[i] = 4;
    162145            goto next_source;
    163146        }           
    164147
    165         test->metric->data.F64[i] = obsMag - fitMag;
    166         test->fitMag->data.F64[i] = fitMag;
     148        try->metric->data.F64[i] = obsMag - fitMag;
     149        try->fitMag->data.F64[i] = fitMag;
    167150        Npsf ++;
    168151
    … …  
    171154
    172155    }
    173     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    174     psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
     156    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     157    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    175158
    176159    // make this optional
    177     // DumpModelFits (test->modelPSF, "modelsPSF.dat");
    178 
    179     // XXX this function wants aperture radius from pmSourcePhotometry
    180     pmPSFMetricModel (test, RADIUS);
    181     psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    182               modelName, test->ApResid, test->dApResid, test->skyBias);
    183 
    184     return (test);
    185 }
    186 
    187 // input: an array of pmModels, pre-allocated psf
    188 // some of the array entries may be NULL, ignore them
    189 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
    190 
    191     // construct the fit vectors from the collection of objects
    192     // use the mask to ignore missing fits
    193     psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
    194     psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
    195     psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
    196     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
    197 
    198     for (int i = 0; i < models->n; i++) {
    199         pmModel *model = models->data[i];
    200         if (model == NULL) continue;
    201 
    202         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
    203         x->data.F64[i] = model->params->data.F32[2];
    204         y->data.F64[i] = model->params->data.F32[3];
    205     }
    206 
    207     // we are doing a robust fit.  after each pass, we drop points which are
    208     // more deviant than three sigma.
    209     // mask is currently updated for each pass. this is inconsistent?
    210 
    211     for (int i = 0; i < psf->params->n; i++) {
    212         for (int j = 0; j < models->n; j++) {
    213             pmModel *model = models->data[j];
    214             if (model == NULL) continue;
    215             z->data.F64[j] = model->params->data.F32[i + 4];
    216             dz->data.F64[j] = 1;
    217             // XXX EAM : need to use actual errors?
    218             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
    219         }
    220 
    221         psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
    222         // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    223         // psPolynomial2DDump (psf->params->data[i]);
    224     }
    225    
    226     psFree (x);
    227     psFree (y);
    228     psFree (z);
    229     psFree (dz);
    230     return (true);
    231 }
    232 
    233 pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) {
    234 
    235     // need to define the relationship between the modelFLT and modelPSF ?
    236    
    237     // find function used to set the model parameters
    238     pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
    239 
    240     // do we need a different model for the PSF vs FLT version?
    241     pmModel *modelPSF = pmModelAlloc (psf->type);
    242 
    243     // set model parameters for this source based on PSF information
    244     modelFromPSFFunc (modelPSF, modelFLT, psf);
    245 
    246     return (modelPSF);
    247 }
    248 
    249 bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
     160    // DumpModelFits (try->modelPSF, "modelsPSF.dat");
     161
     162    // XXX this function wants aperture radius for pmSourcePhotometry
     163    pmPSFtryMetric (try, RADIUS);
     164    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
     165              modelName, try->ApResid, try->dApResid, try->skyBias);
     166
     167    return (try);
     168}
     169
     170
     171bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
    250172
    251173  float dBin;
    252174  int   nKeep, nSkip;
    253175
    254   // the measured (aperture - fit) magnitudes (dA == test->metric)
     176  // the measured (aperture - fit) magnitudes (dA == try->metric)
    255177  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    256178  //     dA = dAo + dsky/flux
    … …  
    262184
    263185  // rflux = ten(0.4*fitMag);
    264   psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
    265   for (int i = 0; i < test->sources->n; i++) {
    266     if (test->mask->data.U8[i]) continue;
    267     rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
     186  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     187  for (int i = 0; i < try->sources->n; i++) {
     188    if (try->mask->data.U8[i]) continue;
     189    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
    268190  }
    269191
    270192  // find min and max of (1/flux):
    271193  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    272   psVectorStats (stats, rflux, NULL, test->mask, 0xff);
     194  psVectorStats (stats, rflux, NULL, try->mask, 0xff);
    273195 
    274196  // build binned versions of rflux, metric
    … …  
    284206  for (int i = 0; i < daBin->n; i++) {
    285207
    286     psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
     208    psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    287209    tmp->n = 0;
    288210
    289211    // accumulate data within bin range
    290     for (int j = 0; j < test->sources->n; j++) {
     212    for (int j = 0; j < try->sources->n; j++) {
    291213      // masked for: bad model fit, outlier in parameters
    292       if (test->mask->data.U8[j]) continue;
     214      if (try->mask->data.U8[j]) continue;
    293215   
    294216      // skip points with extreme dA values
    295       if (fabs(test->metric->data.F64[j]) > 0.5) continue;
     217      if (fabs(try->metric->data.F64[j]) > 0.5) continue;
    296218
    297219      // skip points outside of this bin
    … …  
    299221      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
    300222
    301       tmp->data.F64[tmp->n] = test->metric->data.F64[j];
     223      tmp->data.F64[tmp->n] = try->metric->data.F64[j];
    302224      tmp->n ++;
    303225    }
    … …  
    342264
    343265  // linear fit to rfBin, daBin
    344   psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     266  psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     267
     268  // XXX EAM : this is the intended API (cycle 7? cycle 8?)
     269  // poly = psVectorFitPolynomial1D (poly, maskB, 1, daBin, NULL, rfBin);
     270
     271  // XXX EAM : replace this when the above version is implemented
    345272  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
    346273
    347   psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
     274  psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    348275  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    349276
    … …  
    351278  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    352279
    353   test->ApResid = poly->coeff[0];
    354   test->dApResid = stats->clippedStdev;
    355   test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     280  try->ApResid = poly->coeff[0];
     281  try->dApResid = stats->clippedStdev;
     282  try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    356283
    357284  psFree (rflux);
Note: See TracChangeset for help on using the changeset viewer.