Changeset 33415 for branches/meh_branches/ppstack_test/psphot/src/psphotModelTest.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotModelTest.c (modified) (1 diff)

branches/meh_branches/ppstack_test/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

branches/meh_branches/ppstack_test/psphot/src/psphotModelTest.c

@@ -1 @@
-# include "psphotInternal.h"
-# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
+# include "psphotStandAlone.h"
 
-// XXX add more test information?
-bool psphotModelTest (pmConfig *config, const pmFPAview *view, const char *filerule, psMetadata *recipe) {
+int main (int argc, char **argv) {
 
-    bool status;
-    int modelType = -1;
-    float obsMag, fitMag, value;
-    char name[64];
-    pmPSF *psf = NULL;
-    pmSourceFitMode fitMode;
+    psMemInit();
+    psTimerStart ("complete");
+    pmErrorRegister();                  // register psModule's error codes/messages
+    psphotInit();
 
-    // bit-masks to test for good/bad pixels
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
+    // load command-line arguments, options, and system config data
+    pmConfig *config = psphotModelTestArguments (argc, argv);
+    assert(config);
 
-    // bit-mask to mark pixels not used in analysis
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
+    psphotVersionPrint();
 
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // run model fitting tests on a single source?
-    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
-
-    psTimerStart ("modelTest");
-
-    // find the currently selected readout
-    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
-    PS_ASSERT_PTR_NON_NULL (readout, false);
-
-    // use poissonian errors or local-sky errors
-    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule);
-    if (!status) POISSON_ERRORS = true;
-    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
-
-    // find the various fitting parameters (try test values first)
-    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
-    if (!status || !isfinite(INNER)) {
-        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    }
-    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
-    if (!status || !isfinite(OUTER)) {
-        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    }
-    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
-    if (!status || !isfinite(RADIUS)) {
-        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    }
-    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
-    if (!status || !isfinite(mRADIUS)) {
-        mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    // load input data (config and images (signal, noise, mask)
+    if (!psphotParseCamera (config)) {
+        psErrorStackPrint(stderr, "Error setting up the camera\n");
+        exit (psphotGetExitStatus());
     }
 
-    // define the source of interest
-    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
-    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
-    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
-
-    // what fitting mode to use?
-    fitMode = PM_SOURCE_FIT_EXT;
-    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
-    if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
-    if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
-    if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
-
-    // construct the source structures
-    pmSource *source = pmSourceAlloc();
-    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
-    pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
-
-    // in fitMode, psf sets the model type
-    if (fitMode == PM_SOURCE_FIT_PSF) {
-        psf = psphotLoadPSF (config, view, recipe);
-        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
-        modelType = psf->type;
-        source->type = PM_SOURCE_TYPE_STAR;
-    }
-    if (fitMode == PM_SOURCE_FIT_EXT) {
-        // find the model: supplied by user or first in the PSF_MODEL list
-        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
-        if (!status || !strcasecmp (modelName, "DEFAULT")) {
-            // get the list pointers for the PSF_MODEL entries
-
-            psList *list = NULL;
-            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
-            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
-            if (mdi->type == PS_DATA_STRING) {
-                list = psListAlloc(NULL);
-                psListAdd (list, PS_LIST_HEAD, mdi);
-            } else {
-                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
-                list = psMemIncrRefCounter(mdi->data.list);
-            }
-
-            // take the first list element
-            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
-            modelName = item->data.V;
-        }
-        modelType = pmModelClassGetType (modelName);
-        if (modelType < 0) psAbort("unknown model %s", modelName);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
-    }
-    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
-        // we need to load BOTH a psf and an ext model
-        psf = psphotLoadPSF (config, view, recipe);
-        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
-
-        // find the model: supplied by user or first in the PSF_MODEL list
-        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
-        if (!status || !strcasecmp (modelName, "DEFAULT")) {
-            // get the list pointers for the PSF_MODEL entries
-
-            psList *list = NULL;
-            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
-            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
-            if (mdi->type == PS_DATA_STRING) {
-                list = psListAlloc(NULL);
-                psListAdd (list, PS_LIST_HEAD, mdi);
-            } else {
-                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
-                list = psMemIncrRefCounter(mdi->data.list);
-            }
-
-            // take the first list element
-            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
-            modelName = item->data.V;
-        }
-        modelType = pmModelClassGetType (modelName);
-        if (modelType < 0) psAbort("unknown model %s", modelName);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
+    // call psphot for each readout
+    if (!psphotImageLoop (config, PSPHOT_MODEL_TEST)) {
+        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
+        exit (psphotGetExitStatus());
     }
 
-    // find the local sky
-    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-    if (!status) psAbort("pmSourceLocalSky error");
+    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
 
-    // get the source moments
-    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0, maskVal);
-    if (!status) psAbort("psSourceMoments error");
-    source->peak->value = source->moments->Peak;
+    psErrorCode exit_status = psphotGetExitStatus();
+    psphotCleanup (config);
+    exit (exit_status);
+}
 
-    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->Mx, source->moments->My);
-    fprintf (stderr, "moments: %f, %f - %f\n", source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
-
-    psEllipseMoments moments;
-    moments.x2 = source->moments->Mxx;
-    moments.y2 = source->moments->Myy;
-    moments.xy = source->moments->Mxy;
-    psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
-
-    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
-
-    // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
-    source->modelEXT = model;
-
-    // if any parameters are defined by the user, take those values
-    int nParams = pmModelClassParameterCount (modelType);
-    psF32 *params = model->params->data.F32;
-    params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
-    params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
-    for (int i = 0; i < nParams; i++) {
-        if (i == PM_PAR_XPOS) continue;
-        if (i == PM_PAR_YPOS) continue;
-
-        sprintf (name, "TEST_FIT_PAR%d", i);
-        value = psMetadataLookupF32 (&status, recipe, name);
-        if (status && isfinite (value)) {
-            params[i] = value;
-        }
-    }
-
-    float area = params[4]*params[5];
-    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
-
-    // for PSF fitting, set the shape parameters based on the PSF & source position
-    if (fitMode == PM_SOURCE_FIT_PSF) {
-        source->modelPSF = pmModelFromPSF (model, psf);
-        psFree (model);
-        model = source->modelPSF;
-        params = model->params->data.F32;
-    }
-
-    // list model input shape
-    psEllipseShape shape;
-    shape.sx  = 1.4 / model->params->data.F32[4];
-    shape.sy  = 1.4 / model->params->data.F32[5];
-    shape.sxy = model->params->data.F32[6];
-    axes = psEllipseShapeToAxes (shape, 20.0);
-
-    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
-
-    fprintf (stderr, "input parameters: \n");
-    for (int i = 0; i < nParams; i++) {
-        fprintf (stderr, "%d : %f\n", i, params[i]);
-    }
-
-    // define the pixels used for the fit
-    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", markVal);
-    psphotSaveImage (NULL, source->maskObj, "mask1.fits");
-
-    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
-    if (status) {
-        status = psphotFitSet (source, model, fitset, fitMode, maskVal);
-        exit (0);
-    }
-
-    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
-        // build the psf for the object
-        source->modelPSF = pmModelFromPSF (model, psf);
-        source->modelEXT = model;
-
-        // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
-        int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
-        assert (status);
-
-        model = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
-        params = model->params->data.F32;
-    } else {
-        status = pmSourceFitModel (source, model, fitMode, maskVal);
-    }
-
-    // measure the source mags
-    pmSourcePhotometryModel (&fitMag, model);
-    pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
-    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
-
-    // write out positive object
-    psphotSaveImage (NULL, source->pixels, "object.fits");
-
-    // subtract object, leave local sky
-    // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
-    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-    fprintf (stderr, "output parameters: \n");
-    for (int i = 0; i < nParams; i++) {
-        fprintf (stderr, "%d : %f\n", i, params[i]);
-    }
-
-    // write out
-    psphotSaveImage (NULL, source->pixels, "resid.fits");
-    psphotSaveImage (NULL, source->maskObj, "mask.fits");
-
-    psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));
-
-    exit (0);
-}
+// all functions which return to this level must raise one of the top-level error codes if they
+// exit with an error.  these error codes are used to specify the program exit status

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest