IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 28, 2006, 10:23:51 AM (20 years ago)
Author:
magnier
Message:

This is an important change to the PSF shape model.

  • I have changed the definition of the elliptical contour parameters

(PM_PAR_SXX,SYY are now sigma_x,y/sqrt(2), rather than
sqrt(2)/sigma_x.

  • I have changed the way in which the SXY term is modelled. rather

than fit SXY to a polynomial, I now fit SXY / (1/SXX2 + 1/SYY2)2.
this representation is more likely to be well-described by a
polynomial. the other form is likely to vary as xy/rN rather than xy
rN.

these changes imply that the models all need to treat the SXX,SYY,SXY
terms in the same fashion (ie, fitting a contour of the form (x/SXX)2
+ (y/SYY)2 + SXY*x*y). If an arbitrary model uses some other meaning
for SXX, SYY, SXY, then it will have to define its own exceptions in
the functions pmPSFFromModels and pmModelFromPSF (the latter has an
implementation for each model already).

  • I have cleaned up the pmModel_NAME.c functions to use the #define names

for the parameters (PM_PAR_XXX) everywhere, and to use #def names for
the functions. this latter change makes it easy to specify the model
function names in a single location.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9730 r9770  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-10-24 22:55:05 $
     8 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-10-28 20:23:51 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6969    psFree (psf->ApTrend);
    7070    psFree (psf->growth);
    71     psFree (psf->params);
     71    psFree (psf->params_NEW);
    7272    return;
    7373}
    7474
    75 
    76 
    7775/*****************************************************************************
    78 pmPSFAlloc (type): allocate a pmPSF.
    79     NOTE: a PSF always has 4 parameters fewer than the equivalent model.
    80       This is because those 4 parameters are
     76 pmPSFAlloc (type): allocate a pmPSF.
     77 
     78 NOTE: PSF model parameters which are not modeled on an image are set to NULL in psf->params.
     79 
     80 These are normally:
     81 
    8182 X-center
    8283 Y-center
    … …  
    119120        return(NULL);
    120121    }
    121 
    122     psf->params = psArrayAlloc(Nparams - 4);
    123 
    124     // the order of the PSF parameter (X,Y) fits is determined by the
    125     // psfTrendMask polynomial (user-specified as in the recipe). the
    126     // masks of psfTrendMask are applied to each parameter.
    127     // if psfTrendMask is NULL, these polynomials are not allocated.
    128     // in this case, the user must set them by hand (as in pmPSFfromMD)
    129     // XXX should we drop the hard-wired '4' above and use NULL to identify
    130     // the parameters which are not fitted.  these could be selected by
    131     // testing for the value of PM_PAR_XPOS, etc.
     122    psf->params_NEW = psArrayAlloc(Nparams);
     123
     124    // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
     125    // (user-specified as in the recipe). the masks of psfTrendMask are applied to each
     126    // parameter.  if psfTrendMask is NULL, these polynomials are not allocated.  in this case,
     127    // the user must set them by hand (as in pmPSFfromMD).  the parameters which are not fitted
     128    // are left as NULL.  these are selected by testing for them by the named values (
     129    // PM_PAR_XPOS, etc)
     130
    132131    if (psfTrendMask) {
    133         for (int i = 0; i < psf->params->n; i++) {
     132        for (int i = 0; i < psf->params_NEW->n; i++) {
     133            if (i == PM_PAR_SKY)
     134                continue;
     135            if (i == PM_PAR_I0)
     136                continue;
     137            if (i == PM_PAR_XPOS)
     138                continue;
     139            if (i == PM_PAR_YPOS)
     140                continue;
     141
    134142            psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY);
    135143            for (int nx = 0; nx < param->nX + 1; nx++) {
    … …  
    138146                }
    139147            }
    140             psf->params->data[i] = param;
     148            psf->params_NEW->data[i] = param;
    141149        }
    142150    }
    … …  
    166174    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
    167175
     176    // construct the x,y terms
    168177    for (int i = 0; i < models->n; i++) {
    169178        pmModel *model = models->data[i];
    … …  
    171180            continue;
    172181
    173         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
    174         // note that the data is provided in the F64 array
    175         x->data.F64[i] = model->params->data.F32[2];
    176         y->data.F64[i] = model->params->data.F32[3];
     182        // use F64 for polynomial fitting
     183        x->data.F64[i] = model->params->data.F32[PM_PAR_XPOS];
     184        y->data.F64[i] = model->params->data.F32[PM_PAR_YPOS];
    177185    }
    178186
    … …  
    183191    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    184192
    185     for (int i = 0; i < psf->params->n; i++) {
     193    // skip the unfitted parameters (X, Y, Io, Sky)
     194    for (int i = 0; i < psf->params_NEW->n; i++) {
     195        if (i == PM_PAR_SKY)
     196            continue;
     197        if (i == PM_PAR_I0)
     198            continue;
     199        if (i == PM_PAR_XPOS)
     200            continue;
     201        if (i == PM_PAR_YPOS)
     202            continue;
     203
     204        // select the per-object fitted data for this parameter
    186205        for (int j = 0; j < models->n; j++) {
    187206            pmModel *model = models->data[j];
    188207            if (model == NULL)
    189208                continue;
    190             z->data.F64[j] = model->params->data.F32[i + 4];
    191             dz->data.F64[j] = 1;
    192             // XXX EAM : need to use actual errors?
    193             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
     209            z->data.F64[j] = model->params->data.F32[i];
     210            dz->data.F64[j] = 1; // use the model-fitted error? or S/N?
     211
     212            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
     213            if (i == PM_PAR_SXY) {
     214                z->data.F64[j] = pmPSF_SXYfromModel (model->params->data.F32);
     215            }
    194216        }
    195 
    196         // XXX EAM : this is the expected API (cycle 7? cycle 8?)
    197         psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y);
    198 
    199         // XXX EAM : drop this when above is implemented...
    200         // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
    201 
     217        psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, mask, 0xff, z, dz, x, y);
    202218        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    203         // psPolynomial2DDump (psf->params->data[i]);
     219
     220        // XXX Test output
     221        psPolynomial2D *poly = psf->params_NEW->data[i];
     222        fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
     223        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][0], poly->coeff[1][0], poly->coeff[0][1]);
     224        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][2], poly->coeff[1][1], poly->coeff[2][0]);
     225    }
     226
     227    // XXX test dump of star parameters vs position (compare with fitted values)
     228    {
     229        FILE *f = fopen ("params.dat", "w");
     230
     231        for (int j = 0; j < models->n; j++)
     232        {
     233            pmModel *model = models->data[j];
     234            if (model == NULL)
     235                continue;
     236
     237            pmModel *modelPSF = pmModelFromPSF (model, psf);
     238
     239            fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     240
     241            for (int i = 0; i < psf->params_NEW->n; i++) {
     242                if (psf->params_NEW->data[i] == NULL)
     243                    continue;
     244                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
     245            }
     246            fprintf (f, "%f %d\n", model->chisq, model->nIter);
     247        }
     248        fclose (f);
    204249    }
    205250
    … …  
    211256    return (true);
    212257}
    213 
    214 
    215258
    216259/*****************************************************************************
    … …  
    250293    }
    251294    return true;
     295}
     296
     297// the PSF models the \sigma_{xy} variation of the elliptical contour as a function of position in the image with a
     298// polynomial.  an individual object has a contour of the form (x^2/2sx^2) + (y^2/2sy^2) + sxy*x*y
     299// these are the values of the model->params.  the psf->params term for sxy is actually fitted
     300// to sxy/(sxx^-2 + syy^-2)^2
     301
     302// input: model->param, output: psf->param[PM_PAR_SXY]
     303double pmPSF_SXYfromModel (psF32 *modelPar)
     304{
     305
     306    double SXX = modelPar[PM_PAR_SXX];
     307    double SYY = modelPar[PM_PAR_SYY];
     308    double SXY = modelPar[PM_PAR_SXY];
     309
     310    double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
     311    return (par);
     312}
     313
     314// input: fitted psf->param, output: model->param[PM_PAR_SXY]
     315double pmPSF_SXYtoModel (psF32 *fittedPar)
     316{
     317
     318    double SXX = fittedPar[PM_PAR_SXX];
     319    double SYY = fittedPar[PM_PAR_SYY];
     320    double fit = fittedPar[PM_PAR_SXY];
     321
     322    double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
     323    return SXY;
    252324}
    253325
Note: See TracChangeset for help on using the changeset viewer.