IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32154


Ignore:
Timestamp:
Aug 21, 2011, 10:30:23 AM (15 years ago)
Author:
eugene
Message:

fixed up the psphotModelTest program (fit a single source)

Location:
branches/eam_branches/ipp-20110710/psphot/src
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110710/psphot/src/Makefile.am

    r31989 r32154  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack
     27bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack psphotModelTest
    2828# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy
    2929
    … …  
    4343psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    4444psphotStack_LDADD = libpsphot.la
     45
     46psphotModelTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     47psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     48psphotModelTest_LDADD = libpsphot.la
    4549
    4650# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    100104        psphotStackObjects.c          \
    101105        psphotStackPSF.c              \
     106        psphotCleanup.c
     107
     108# a psphot-variant that simply generates the PSF model
     109psphotModelTest_SOURCES = \
     110        psphotModelTest.c            \
     111        psphotModelTestArguments.c   \
     112        psphotParseCamera.c        \
     113        psphotImageLoop.c   \
     114        psphotMosaicChip.c         \
    102115        psphotCleanup.c
    103116
    … …  
    129142        psphotMakePSFReadout.c         \
    130143        psphotModelBackground.c        \
     144        psphotModelTestReadout.c       \
    131145        psphotMaskBackground.c         \
    132146        psphotSubtractBackground.c     \
    … …  
    192206        psphotEfficiency.c
    193207
    194 # XXX need to fix this for the new apis
    195 #       psphotModelTest.c             
    196 
    197208# re-instate these
    198209#       psphotIsophotal.c              \

  • branches/eam_branches/ipp-20110710/psphot/src/psphotModelTest.c

    r31452 r32154  
    1 # include "psphotInternal.h"
    2 # define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
     1# include "psphotStandAlone.h"
    32
    4 // XXX add more test information?
    5 bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) {
     3int main (int argc, char **argv) {
    64
    7     bool status;
    8     int modelType = -1;
    9     float obsMag, fitMag, value;
    10     char name[64];
    11     pmPSF *psf = NULL;
    12     pmSourceFitMode fitMode;
     5    psMemInit();
     6    psTimerStart ("complete");
     7    pmErrorRegister();                  // register psModule's error codes/messages
     8    psphotInit();
    139
    14     // bit-masks to test for good/bad pixels
    15     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    16     assert (maskVal);
     10    // load command-line arguments, options, and system config data
     11    pmConfig *config = psphotModelTestArguments (argc, argv);
     12    assert(config);
    1713
    18     // bit-mask to mark pixels not used in analysis
    19     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    20     assert (markVal);
     14    psphotVersionPrint();
    2115
    22     // maskVal is used to test for rejected pixels, and must include markVal
    23     maskVal |= markVal;
    24 
    25     // run model fitting tests on a single source?
    26     if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
    27 
    28     psTimerStart ("modelTest");
    29 
    30     // find the currently selected readout
    31     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    32     PS_ASSERT_PTR_NON_NULL (readout, false);
    33 
    34     // use poissonian errors or local-sky errors
    35     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule);
    36     if (!status) POISSON_ERRORS = true;
    37     pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
    38 
    39     // find the various fitting parameters (try test values first)
    40     float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    41     if (!status || !isfinite(INNER)) {
    42         INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    43     }
    44     float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    45     if (!status || !isfinite(OUTER)) {
    46         OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    47     }
    48     float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    49     if (!status || !isfinite(RADIUS)) {
    50         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    51     }
    52     float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    53     if (!status || !isfinite(mRADIUS)) {
    54         mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     16    // load input data (config and images (signal, noise, mask)
     17    if (!psphotParseCamera (config)) {
     18        psErrorStackPrint(stderr, "Error setting up the camera\n");
     19        exit (psphotGetExitStatus());
    5520    }
    5621
    57     // define the source of interest
    58     float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    59     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
    60     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
    61 
    62     // what fitting mode to use?
    63     fitMode = PM_SOURCE_FIT_EXT;
    64     char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    65     if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
    66     if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
    67     if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
    68 
    69     // construct the source structures
    70     pmSource *source = pmSourceAlloc();
    71     source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
    72     pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
    73 
    74     // in fitMode, psf sets the model type
    75     if (fitMode == PM_SOURCE_FIT_PSF) {
    76         psf = psphotLoadPSF (config, view, recipe);
    77         if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    78         modelType = psf->type;
    79         source->type = PM_SOURCE_TYPE_STAR;
    80     }
    81     if (fitMode == PM_SOURCE_FIT_EXT) {
    82         // find the model: supplied by user or first in the PSF_MODEL list
    83         char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    84         if (!status || !strcasecmp (modelName, "DEFAULT")) {
    85             // get the list pointers for the PSF_MODEL entries
    86 
    87             psList *list = NULL;
    88             psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
    89             if (mdi == NULL) psAbort("missing PSF_MODEL selection");
    90             if (mdi->type == PS_DATA_STRING) {
    91                 list = psListAlloc(NULL);
    92                 psListAdd (list, PS_LIST_HEAD, mdi);
    93             } else {
    94                 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
    95                 list = psMemIncrRefCounter(mdi->data.list);
    96             }
    97 
    98             // take the first list element
    99             psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    100             modelName = item->data.V;
    101         }
    102         modelType = pmModelClassGetType (modelName);
    103         if (modelType < 0) psAbort("unknown model %s", modelName);
    104         source->type = PM_SOURCE_TYPE_EXTENDED;
    105     }
    106     if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
    107         // we need to load BOTH a psf and an ext model
    108         psf = psphotLoadPSF (config, view, recipe);
    109         if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    110 
    111         // find the model: supplied by user or first in the PSF_MODEL list
    112         char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    113         if (!status || !strcasecmp (modelName, "DEFAULT")) {
    114             // get the list pointers for the PSF_MODEL entries
    115 
    116             psList *list = NULL;
    117             psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
    118             if (mdi == NULL) psAbort("missing PSF_MODEL selection");
    119             if (mdi->type == PS_DATA_STRING) {
    120                 list = psListAlloc(NULL);
    121                 psListAdd (list, PS_LIST_HEAD, mdi);
    122             } else {
    123                 if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
    124                 list = psMemIncrRefCounter(mdi->data.list);
    125             }
    126 
    127             // take the first list element
    128             psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    129             modelName = item->data.V;
    130         }
    131         modelType = pmModelClassGetType (modelName);
    132         if (modelType < 0) psAbort("unknown model %s", modelName);
    133         source->type = PM_SOURCE_TYPE_EXTENDED;
     22    // call psphot for each readout
     23    if (!psphotImageLoop (config, PSPHOT_MODEL_TEST)) {
     24        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
     25        exit (psphotGetExitStatus());
    13426    }
    13527
    136     // find the local sky
    137     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    138     if (!status) psAbort("pmSourceLocalSky error");
     28    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
    13929
    140     // get the source moments
    141     status = pmSourceMoments (source, mRADIUS, 0.0, 1.0, maskVal);
    142     if (!status) psAbort("psSourceMoments error");
    143     source->peak->value = source->moments->Peak;
     30    psErrorCode exit_status = psphotGetExitStatus();
     31    psphotCleanup (config);
     32    exit (exit_status);
     33}
    14434
    145     fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
    146     fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
    147 
    148     psEllipseMoments moments;
    149     moments.x2 = source->moments->Mxx;
    150     moments.y2 = source->moments->Myy;
    151     moments.xy = source->moments->Mxy;
    152     psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
    153 
    154     fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    155 
    156     // get the initial model parameter guess
    157     pmModel *model = pmSourceModelGuess (source, modelType);
    158     source->modelEXT = model;
    159 
    160     // if any parameters are defined by the user, take those values
    161     int nParams = pmModelClassParameterCount (modelType);
    162     psF32 *params = model->params->data.F32;
    163     params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
    164     params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
    165     for (int i = 0; i < nParams; i++) {
    166         if (i == PM_PAR_XPOS) continue;
    167         if (i == PM_PAR_YPOS) continue;
    168 
    169         sprintf (name, "TEST_FIT_PAR%d", i);
    170         value = psMetadataLookupF32 (&status, recipe, name);
    171         if (status && isfinite (value)) {
    172             params[i] = value;
    173         }
    174     }
    175 
    176     float area = params[4]*params[5];
    177     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
    178 
    179     // for PSF fitting, set the shape parameters based on the PSF & source position
    180     if (fitMode == PM_SOURCE_FIT_PSF) {
    181         source->modelPSF = pmModelFromPSF (model, psf);
    182         psFree (model);
    183         model = source->modelPSF;
    184         params = model->params->data.F32;
    185     }
    186 
    187     // list model input shape
    188     psEllipseShape shape;
    189     shape.sx  = 1.4 / model->params->data.F32[4];
    190     shape.sy  = 1.4 / model->params->data.F32[5];
    191     shape.sxy = model->params->data.F32[6];
    192     axes = psEllipseShapeToAxes (shape, 20.0);
    193 
    194     fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    195 
    196     fprintf (stderr, "input parameters: \n");
    197     for (int i = 0; i < nParams; i++) {
    198         fprintf (stderr, "%d : %f\n", i, params[i]);
    199     }
    200 
    201     // define the pixels used for the fit
    202     psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", markVal);
    203     psphotSaveImage (NULL, source->maskObj, "mask1.fits");
    204 
    205     char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    206     if (status) {
    207         status = psphotFitSet (source, model, fitset, fitMode, maskVal);
    208         exit (0);
    209     }
    210 
    211     if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
    212         // build the psf for the object
    213         source->modelPSF = pmModelFromPSF (model, psf);
    214         source->modelEXT = model;
    215 
    216         // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
    217         int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    218         assert (status);
    219 
    220         model = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
    221         params = model->params->data.F32;
    222     } else {
    223         status = pmSourceFitModel (source, model, fitMode, maskVal);
    224     }
    225 
    226     // measure the source mags
    227     pmSourcePhotometryModel (&fitMag, model);
    228     pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
    229     fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
    230 
    231     // write out positive object
    232     psphotSaveImage (NULL, source->pixels, "object.fits");
    233 
    234     // subtract object, leave local sky
    235     // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
    236     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    237 
    238     fprintf (stderr, "output parameters: \n");
    239     for (int i = 0; i < nParams; i++) {
    240         fprintf (stderr, "%d : %f\n", i, params[i]);
    241     }
    242 
    243     // write out
    244     psphotSaveImage (NULL, source->pixels, "resid.fits");
    245     psphotSaveImage (NULL, source->maskObj, "mask.fits");
    246 
    247     psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));
    248 
    249     exit (0);
    250 }
     35// all functions which return to this level must raise one of the top-level error codes if they
     36// exit with an error.  these error codes are used to specify the program exit status
Note: See TracChangeset for help on using the changeset viewer.