IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:32:31 PM (15 years ago)
Author:
eugene
Message:
  • clarified the output to show the stages and their timing more clearly
  • psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
  • psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
  • for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
  • extended source analysis (petrosians) is threaded
  • iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
  • the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
  • dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
  • fake sources inserted for the efficiency analysis are removed at the end of the analysis
  • radial apertures is now threaded
  • edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
  • new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
  • define the fit region and window region based on the 1st radial moments
  • use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
  • remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis
Location:
trunk/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        2222psphotMakePSF
        2323psphotStack
         24psphotModelTest

  • trunk/psphot/src/psphotModelTest.c

    r31452 r32348  
    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.