IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 20, 2007, 10:34:30 AM (19 years ago)
Author:
eugene
Message:

finished PSFConvModel code; added structure to carry convolved model data

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotModelTest.c

    r14326 r14338  
    11# include "psphotInternal.h"
    2 static char DEFAULT_MODE[] = "EXT";
    3 
    4 // XXX consider this function : add more test information?
    5 bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
     2# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
     3
     4// XXX add more test information?
     5bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
    66
    77    bool status;
    88    int modelType;
    9     unsigned int Nfail;
    109    float obsMag, fitMag, value;
    1110    char name[64];
    … …  
    1312    pmSourceFitMode fitMode;
    1413
    15     psMetadataItem *item  = NULL;
     14    // run model fitting tests on a single source?
     15    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
     16
     17    // find the currently selected readout
     18    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
     19    PS_ASSERT_PTR_NON_NULL (readout, false);
    1620
    1721    // use poissonian errors or local-sky errors
    … …  
    2024    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
    2125
    22     // run model fitting tests on a single source
    23     if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
     26    // find the various fitting parameters (try test values first)
     27    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
     28    if (!status || !isfinite(INNER)) {
     29        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     30    }
     31    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
     32    if (!status || !isfinite(OUTER)) {
     33        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     34    }
     35    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
     36    if (!status || !isfinite(RADIUS)) {
     37        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
     38    }
     39    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
     40    if (!status || !isfinite(mRADIUS)) {
     41        mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     42    }
     43
     44    // define the source of interest
     45    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     46    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     47    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
    2448
    2549    // what fitting mode to use?
     50    fitMode = PM_SOURCE_FIT_EXT;
    2651    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    27     if (!status || !strcasecmp (fitModeWord, "DEFAULT")) {
    28         fitModeWord = DEFAULT_MODE;
    29     }
    30     fitMode = PM_SOURCE_FIT_EXT;
    31     if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
     52    if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
     53    if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
     54    if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
    3255
    3356    // in fitMode, psf sets the model type
    3457    if (fitMode == PM_SOURCE_FIT_PSF) {
    35         // XXX load psf using psphotLoadPSF
    36         char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
    37         if (!status) psAbort("PSF_INPUT_FILE not supplied");
    38         psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE);
    39         psf = pmPSFfromMetadata (psfData);
     58        psf = psphotLoadPSF (config, view, recipe);
     59        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    4060        modelType = psf->type;
    41     } else {
     61    }
     62    if (fitMode == PM_SOURCE_FIT_EXT) {
    4263        // find the model: supplied by user or first in the PSF_MODEL list
    4364        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    … …  
    5778
    5879            // take the first list element
    59             item = psListGet (list, PS_LIST_HEAD);
     80            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    6081            modelName = item->data.V;
    6182        }
    … …  
    6384        if (modelType < 0) psAbort("unknown model %s", modelName);
    6485    }
    65 
    66     // find the fitting parameters (try test values first)
    67     float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    68     if (!status || !isfinite(INNER)) {
    69         INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    70     }
    71 
    72     float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    73     if (!status || !isfinite(OUTER)) {
    74         OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    75     }
    76 
    77     float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    78     if (!status || !isfinite(RADIUS)) {
    79         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    80     }
    81 
    82     float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    83     if (!status || !isfinite(mRADIUS)) {
    84         mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    85     }
    86 
    87     // define the source of interest
    88     float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    89     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     86    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
     87        // we need to load BOTH a psf and an ext model
     88        psf = psphotLoadPSF (config, view, recipe);
     89        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
     90
     91        // find the model: supplied by user or first in the PSF_MODEL list
     92        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
     93        if (!status || !strcasecmp (modelName, "DEFAULT")) {
     94            // get the list pointers for the PSF_MODEL entries
     95
     96            psList *list = NULL;
     97            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
     98            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
     99            if (mdi->type == PS_DATA_STRING) {
     100                list = psListAlloc(NULL);
     101                psListAdd (list, PS_LIST_HEAD, mdi);
     102            } else {
     103                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
     104                list = psMemIncrRefCounter(mdi->data.list);
     105            }
     106
     107            // take the first list element
     108            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
     109            modelName = item->data.V;
     110        }
     111        modelType = pmModelSetType (modelName);
     112        if (modelType < 0) psAbort("unknown model %s", modelName);
     113    }
    90114
    91115    // construct the source structures
    … …  
    100124    // get the source moments
    101125    status = pmSourceMoments (source, mRADIUS);
    102     if (!status) psAbort("pmSourceLocalSky error");
     126    if (!status) psAbort("psSourceMoments error");
    103127    source->peak->value = source->moments->Peak;
    104128
    … …  
    116140    // get the initial model parameter guess
    117141    pmModel *model = pmSourceModelGuess (source, modelType);
    118     // if any parameters are defined, use those values
     142
     143    // if any parameters are defined by the user, take those values
    119144    int nParams = pmModelParameterCount (modelType);
    120145    psF32 *params = model->params->data.F32;
     146    params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
     147    params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
    121148    for (int i = 0; i < nParams; i++) {
    122         if (i == 2) {
    123             params[i] = xObj;
    124             continue;
    125         }
    126         if (i == 3) {
    127             params[i] = yObj;
    128             continue;
    129         }
     149        if (i == PM_PAR_XPOS) continue;
     150        if (i == PM_PAR_YPOS) continue;
     151
    130152        sprintf (name, "TEST_FIT_PAR%d", i);
    131153        value = psMetadataLookupF32 (&status, recipe, name);
    … …  
    138160    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
    139161
     162    // for PSF fitting, set the shape parameters based on the PSF & source position
    140163    if (fitMode == PM_SOURCE_FIT_PSF) {
    141164        pmModel *modelPSF = pmModelFromPSF (model, psf);
    … …  
    154177    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    155178
    156 
    157179    fprintf (stderr, "input parameters: \n");
    158180    for (int i = 0; i < nParams; i++) {
    … …  
    162184    // define the pixels used for the fit
    163185    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", mark);
     186    psphotSaveImage (NULL, source->maskObj, "mask1.fits");
    164187
    165188    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    … …  
    169192    }
    170193
    171     status = pmSourceFitModel (source, model, fitMode, maskVal);
     194    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
     195        // build the psf for the object
     196        source->modelPSF = pmModelFromPSF (model, psf);
     197        source->modelEXT = model;
     198        status = psphotPSFConvModel (source, recipe, maskVal);
     199        model = source->modelConv;
     200        params = model->params->data.F32;
     201    } else {
     202        status = pmSourceFitModel (source, model, fitMode, maskVal);
     203    }
    172204
    173205    // measure the source mags
Note: See TracChangeset for help on using the changeset viewer.