Index: /branches/eam_branches/ipp-20110710/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/Makefile.am (revision 32153) +++ /branches/eam_branches/ipp-20110710/psphot/src/Makefile.am (revision 32154) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack +bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack psphotModelTest # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -43,4 +43,8 @@ psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotStack_LDADD = libpsphot.la + +psphotModelTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotModelTest_LDADD = libpsphot.la # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -100,4 +104,13 @@ psphotStackObjects.c \ psphotStackPSF.c \ + psphotCleanup.c + +# a psphot-variant that simply generates the PSF model +psphotModelTest_SOURCES = \ + psphotModelTest.c \ + psphotModelTestArguments.c \ + psphotParseCamera.c \ + psphotImageLoop.c \ + psphotMosaicChip.c \ psphotCleanup.c @@ -129,4 +142,5 @@ psphotMakePSFReadout.c \ psphotModelBackground.c \ + psphotModelTestReadout.c \ psphotMaskBackground.c \ psphotSubtractBackground.c \ @@ -192,7 +206,4 @@ psphotEfficiency.c -# XXX need to fix this for the new apis -# psphotModelTest.c - # re-instate these # psphotIsophotal.c \ Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTest.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTest.c (revision 32153) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTest.c (revision 32154) @@ -1,250 +1,36 @@ -# 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 Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestArguments.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestArguments.c (revision 32154) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestArguments.c (revision 32154) @@ -0,0 +1,189 @@ +# include "psphotStandAlone.h" + +static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile) +{ + fprintf(ofile, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "\n" + "%s -help\n" + "\n" + "%s -version\n" + "\n" + "where:\n" + " FileNameList is a text file containing filenames, one per line\n" + " MaskFileNameList is a text file of mask filenames, one per line\n" + " VarFileNameList is a text file of variance filenames, one per line\n" + " OutFileBaseName is the 'root name' for output files\n" + "\n" + "additional options:\n" + " -psf PsfFile1[,PsfFile2,...] or -psflist PsfFileNameList\n" + " specify PSF rather than letting %s estimate it\n" + " -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n" + " specify additional sources for PSF generation\n" + " -chip nn[,nn,...]\n" + " select detector chips to process; default is all.\n" + " Indices correspond to zero-based offset in the FPA metadata table.\n" + " -photcode PhotoCodeName\n" + " specify photocode\n" + " -region RegionString\n" + " specify analysis region. String is of form '[x0:x1,y0:y1]'\n" + " To use this option you must define a default in psphot.config\n" + " -visual\n" + " turns on interactive display mode\n" + " -dumpconfig CfgFileName\n" + " causes config info to be dumped to the named file.\n" + " -break NOTHING|BACKMDL|PEAKS|MOMENTS|PSFMODEL|ENSEMBLE|PASS1\n" + " choose a point at which to exit processing early\n" + " -nthreads n\n" + " set number of parallel threads of execution\n" + " -F OldFileRule ReplacementFileRule\n" + " change file naming rule; e.g. '-F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF'\n" + " -D name stringval\n" + " set a string-valued config parameter\n" + " -Di name intval\n" + " set an integer-valued config parameter\n" + " -Df name fval\n" + " set a float-valued config parameter\n" + " -Db name boolval\n" + " set a boolean-valued config parameter\n" + " -v, -vv, -vvv\n" + " set increasing levels of verbosity\n" + " -logfmt FormatString\n" + " set format string used for log messages\n" + " -trace Fac Lvl\n" + " set tracing for facility Fac to integer Lvl, e.g. '-trace err 10'\n" + " -trace-levels\n" + " print current trace levels\n", + program,program,program,program,program); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_SUCCESS); +} + +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config, // Configuration + int exitCode + ) +{ + fprintf(stderr, + "Usage: one of the following\n" + "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n" + " -variance varfile1[,varfile2,...] OutFileBaseName\n" + "\n" + "%s -list FileNameList [-masklist MaskFileNameList] \n" + " -variancelist VarFileNameList OutFileBaseName\n" + "\n" + "Try '%s -help' for more options and explanation\n", + program,program,program); + if (exitCode != PS_EXIT_SUCCESS) + psErrorStackPrint(stderr, "Error reading arguments\n"); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(exitCode); +} + +pmConfig *psphotModelTestArguments(int argc, char **argv) { + + int N; + bool status; + + // load config data from default locations + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + + PS_ARGUMENTS_GENERIC( psphot, config, argc, argv ); + + // save the following additional recipe values based on command-line options + // these options override the PSPHOT recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); + + // Number of threads is handled + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ) + + // visual : interactive display mode + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + pmVisualSetVisual(true); + } + + // break : used from recipe throughout psphotReadout + if ((N = psArgumentGet (argc, argv, "-break"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF + if ((N = psArgumentGet (argc, argv, "-region"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-chip"))) { + if (argc<=N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + exit(PS_EXIT_CONFIG_ERROR); + } + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // if these command-line options are supplied, load the file name lists into config->arguments + // override any configuration-specified source for these files + // + pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); + + // an input list of sources to use is allowed, but not required + pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist"); + pmConfigFileSetsMD (config->arguments, &argc, argv, "SRCTEXT", "-srctext", "-srctextlist"); + + if (argc == 1) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (psArgumentGet(argc, argv, "-help") || + psArgumentGet(argc, argv, "-h")) + writeHelpInfo(argv[0], config, stdout); + + // the input file is a required argument; if not found, we will exit + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); + if (!status) { + psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments"); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + if (argc != 2) { + psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1); + usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR); + } + + // output position is fixed + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psTrace("psphot", 1, "Done with psphotModelTestArguments...\n"); + return (config); +} Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c (revision 32154) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c (revision 32154) @@ -0,0 +1,204 @@ +# include "psphotInternal.h" +# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY + +bool psphotModelTestReadout (pmConfig *config, const pmFPAview *view, const char *filerule) { + + bool status; + + pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // remove cruft from the input analysis structure + if (!psphotCleanInputs (config, view, filerule)) { + psError (PSPHOT_ERR_PROG, false, "trouble setting up the inputs"); + return false; + } + + // set the photcode for this image + if (!psphotAddPhotcode (config, view, filerule)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } + + // Generate the mask and weight images, including the user-defined analysis region of interest + if (!psphotSetMaskAndVariance (config, view, filerule)) { + return psphotReadoutCleanup(config, view, filerule); + } + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + + if (!psphotSubtractBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + + float MIN_KRON_RADIUS = 5.0; + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // use poissonian errors or local-sky errors + bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, filerule); + if (!status) POISSON_ERRORS = true; + + // 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"); + } + + // 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"); + + // construct the source structures + pmSource *source = pmSourceAlloc(); + source->peak = pmPeakAlloc (xObj, yObj, 0, 0); + pmSourceDefinePixels (source, readout, xObj, yObj, OUTER); + + // find the model: supplied by user or first in the PSF_MODEL list + char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); + int modelType = pmModelClassGetType (modelName); + if (modelType < 0) psAbort("unknown model %s", modelName); + source->type = PM_SOURCE_TYPE_EXTENDED; + + // find the local sky + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); + if (!status) psAbort("pmSourceLocalSky error"); + + // get the source moments + status = pmSourceMoments (source, mRADIUS, 0.0, 0.0, MIN_KRON_RADIUS, maskVal); + if (!status) psAbort("psSourceMoments error"); + + 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; + + char name[32]; + sprintf (name, "TEST_FIT_PAR%d", i); + float 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); + + // 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"); + + float SKY_SIG = psMetadataLookupF32(&status, readout->analysis, "SKY_STDEV"); + + // options which modify the behavior of the model fitting + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->nIter = psMetadataLookupS32(&status, recipe, "PSF_FIT_ITER"); // Maximum number of fit iterations + fitOptions->minTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MIN_TOL"); // Fit tolerance + fitOptions->maxTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_TOL"); // Fit tolerance + fitOptions->maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance + fitOptions->poissonErrors = POISSON_ERRORS; + fitOptions->weight = PS_SQR(SKY_SIG); + fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + + status = pmSourceFitModel (source, model, fitOptions, maskVal); + + // measure the source mags + float fitMag = NAN; + float fitFlux = NAN; + float obsMag = NAN; + pmSourcePhotometryModel (&fitMag, &fitFlux, 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); +}