Index: /branches/rel10_ifa/psModules/src/objects/Makefile.am =================================================================== --- /branches/rel10_ifa/psModules/src/objects/Makefile.am (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/Makefile.am (revision 6556) @@ -4,10 +4,24 @@ libpsmoduleobjects_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmoduleobjects_la_SOURCES = \ - pmObjects.c \ + pmPeaks.c \ + pmMoments.c \ + pmModel.c \ + pmModelGroup.c \ + pmSource.c \ + pmSourceSky.c \ + pmSourceContour.c \ + pmSourceFitModel.c \ + pmSourceFitSet.c \ + pmSourcePhotometry.c \ + pmSourceIO.c \ + pmSourceIO_CMF.c \ + pmSourceIO_CMP.c \ + pmSourceIO_OBJ.c \ + pmSourceIO_SX.c \ pmPSF.c \ pmPSFtry.c \ - pmModelGroup.c \ pmGrowthCurve.c \ - psEllipse.c + pmFPAviewReadObjects.c \ + pmFPAviewWriteObjects.c EXTRA_DIST = \ @@ -19,8 +33,16 @@ psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ - pmObjects.h \ + pmPeaks.h \ + pmMoments.h \ + pmModel.h \ + pmModelGroup.h \ + pmSource.h \ + pmSourceSky.h \ + pmSourceContour.h \ + pmSourceFitModel.h \ + pmSourceFitSet.h \ + pmSourcePhotometry.h \ + pmSourceIO.h \ pmPSF.h \ pmPSFtry.h \ - pmModelGroup.h \ - pmGrowthCurve.h \ - psEllipse.h + pmGrowthCurve.h Index: /branches/rel10_ifa/psModules/src/objects/pmFPAviewObjectsIO.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmFPAviewObjectsIO.h (revision 6556) +++ /branches/rel10_ifa/psModules/src/objects/pmFPAviewObjectsIO.h (revision 6556) @@ -0,0 +1,27 @@ +/** @file pmFPAviewObjectsIO.h + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_FPA_OBJECTS_IO_H +# define PM_FPA_OBJECTS_IO_H + +bool pmFPAviewReadObjects (pmFPAview *view, pmFPAfile *file); +bool pmFPAReadObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file); +bool pmChipReadObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file); +bool pmCellReadObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file); +bool pmReadoutReadObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file); + +bool pmFPAviewWriteObjects (pmFPAview *view, pmFPAfile *file); +bool pmFPAWriteObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file); +bool pmChipWriteObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file); +bool pmCellWriteObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file); +bool pmReadoutWriteObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file); + +# endif /* PM_FPA_OBJECTS_IO_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c (revision 6556) @@ -1,3 +1,29 @@ -# include "pmFPAviewReadObjects.h" +/** @file pmSourcePhotometry.c + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" + +#include "pmSourceIO.h" +#include "pmFPAviewObjectsIO.h" // Given a FITS file pointer, read the table of object data @@ -42,10 +68,10 @@ // read in all chip-level Objects files for this FPA -pmFPAReadObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) +bool pmFPAReadObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->data[i]; + pmChip *chip = fpa->chips->data[i]; pmChipReadObjects (chip, view, file); } @@ -54,10 +80,10 @@ // read in all cell-level Objects files for this chip -pmChipReadObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) +bool pmChipReadObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->data[i]; + pmCell *cell = chip->cells->data[i]; pmCellReadObjects (cell, view, file); } @@ -66,10 +92,10 @@ // read in all readout-level Objects files for this cell -pmCellReadObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) +bool pmCellReadObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < cell->readouts->n; i++) { - pmReadout *readout = cell->data[i]; + pmReadout *readout = cell->readouts->data[i]; pmReadoutReadObjects (readout, view, file); } @@ -78,8 +104,10 @@ // read in all readout-level Objects files for this cell -pmReadoutReadObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) +bool pmReadoutReadObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) { + bool status; psArray *sources; + pmHDU *hdu; switch (file->type) { @@ -96,9 +124,9 @@ // read in header, if not yet loaded - pmHDU *hdu = pmFPAviewThisHDU (view); + hdu = pmFPAviewThisHDU (view); - filename = pmConfigNameFromRule (file->filerule, view); + char *filename = pmFPAviewNameFromRule (file->filerule, view); file->fits = psFitsOpen (filename, "r"); - hdu->header = psFitsReadHeader (file->fits); + hdu->header = psFitsReadHeader (NULL, file->fits); psFitsClose (file->fits); sources = pmSourcesReadCMP (filename, hdu->header); @@ -107,14 +135,18 @@ case PM_FPA_FILE_CMF: // read in header, if not yet loaded - pmHDU *hdu = pmFPAviewThisHDU (view); + hdu = pmFPAviewThisHDU (view); if (hdu->header == NULL) { - headname = pmConfigNameFromRule (file->extxtra, view); - psFitsMoveExtName (fits, headname); - hdu->header = psFitsReadHeader (file->fits); + char *headname = pmFPAviewNameFromRule (file->extxtra, view); + psFitsMoveExtName (file->fits, headname); + hdu->header = psFitsReadHeader (NULL, file->fits); } - dataname = pmConfigNameFromRule (file->extrule, view); - psFitsMoveExtName (fits, dataname); + char *dataname = pmFPAviewNameFromRule (file->extrule, view); + psFitsMoveExtName (file->fits, dataname); sources = pmSourcesReadCMF (file->fits, hdu->header); + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); break; } Index: /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c (revision 6556) @@ -1,3 +1,29 @@ -# include "pmFPAviewWriteObjects.h" +/** @file pmSourcePhotometry.c + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" + +#include "pmSourceIO.h" +#include "pmFPAviewObjectsIO.h" // Given a FITS file pointer, read the table of object data @@ -42,10 +68,10 @@ // read in all chip-level Objects files for this FPA -pmFPAWriteObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) +bool pmFPAWriteObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->data[i]; + pmChip *chip = fpa->chips->data[i]; pmChipWriteObjects (chip, view, file); } @@ -54,10 +80,10 @@ // read in all cell-level Objects files for this chip -pmChipWriteObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) +bool pmChipWriteObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < chip->cells->n; i++) { - pmCell *cell = chip->data[i]; + pmCell *cell = chip->cells->data[i]; pmCellWriteObjects (cell, view, file); } @@ -66,10 +92,10 @@ // read in all readout-level Objects files for this cell -pmCellWriteObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) +bool pmCellWriteObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) { for (int i = 0; i < cell->readouts->n; i++) { - pmReadout *readout = cell->data[i]; + pmReadout *readout = cell->readouts->data[i]; pmReadoutWriteObjects (readout, view, file); } @@ -78,6 +104,12 @@ // read in all readout-level Objects files for this cell -pmReadoutWriteObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) +bool pmReadoutWriteObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) { + + bool status; + char *filename; + char *dataname; + char *headname; + pmHDU *hdu; psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); @@ -85,23 +117,23 @@ switch (file->type) { case PM_FPA_FILE_OBJ: - filename = pmConfigNameFromRule (file->filerule, view); - pmSourcesWriteOBJ (filename, sources); + filename = pmFPAviewNameFromRule (file->filerule, view); + pmSourcesWriteOBJ (sources, filename); break; case PM_FPA_FILE_SX: - filename = pmConfigNameFromRule (file->filerule, view); - pmSourcesWriteSX (filename, sources); + filename = pmFPAviewNameFromRule (file->filerule, view); + pmSourcesWriteSX (sources, filename); break; case PM_FPA_FILE_CMP: // a SPLIT format : only one header and object table per file - pmHDU *hdu = pmFPAviewThisHDU (view); - filename = pmConfigNameFromRule (file->filerule, view); - pmSourcesWriteCMP (filename, sources, hdu->header); + hdu = pmFPAviewThisHDU (view); + filename = pmFPAviewNameFromRule (file->filerule, view); + pmSourcesWriteCMP (sources, filename, hdu->header); break; case PM_FPA_FILE_CMF: // write header, if not yet written? - pmHDU *hdu = pmFPAviewThisHDU (view); + hdu = pmFPAviewThisHDU (view); // XXX is this needed? is it automatically added? @@ -111,5 +143,5 @@ if (hdu->header) { psMetadataItem *mdi = NULL; - headname = pmConfigNameFromRule (file->extxtra, view); + headname = pmFPAviewNameFromRule (file->extxtra, view); // set NAXIS to 0 (we don't write out the data array) @@ -119,8 +151,15 @@ mdi->type = PS_DATA_S32; - psFitsWriteHeader (file->fits, hdu->header, headname); + if (strcasecmp (headname, "PHU")) { + psMetadataAdd (hdu->header, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", headname); + } + psFitsWriteHeader (hdu->header, file->fits); } - dataname = pmConfigNameFromRule (file->extrule, view); + dataname = pmFPAviewNameFromRule (file->extrule, view); pmSourcesWriteCMF (file->fits, sources, hdu->header, dataname); + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); break; } Index: /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.c (revision 6556) @@ -1,4 +1,17 @@ +/** @file pmGrowthCurve.c + * + * Measure the curve-of-growth for sources + * + * @author EAM, IfA + * + * @version $Revision: 1.1.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + #include -# include "pmGrowthCurve.h" +#include "pmGrowthCurve.h" static void pmGrowthCurveFree (pmGrowthCurve *growth) Index: /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmGrowthCurve.h (revision 6556) @@ -1,2 +1,3 @@ + # ifndef PM_GROWTH_CURVE_H # define PM_GROWTH_CURVE_H @@ -19,3 +20,5 @@ psF32 pmGrowthCurveCorrect (pmGrowthCurve *growth, psF32 radius); +// XXX psVectorBracket,Interpolate should be put in pslib + # endif /* PM_GROWTH_CURVE_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmModel.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModel.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmModel.c (revision 6556) @@ -1,11 +1,11 @@ -/** @file pmObjects.c - * - * This file will ... +/** @file pmModel.c + * + * Functions to define and manipulate object models * * @author GLG, MHPCC - * @author EAM, IfA: significant modifications. - * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,6 +17,15 @@ #include #include "pslib.h" -#include "pmObjects.h" -#include "pmModelGroup.h" +#include "pmModel.h" + +/* The following types and functions need to be known by pmModel.c and are defined in + pmModelGroup.h. The use of pmSource and other types in pmModelGroup.h prevent us from + directly including it here ***/ + +typedef psMinimizeLMChi2Func pmModelFunc; +typedef bool (*pmModelFitStatusFunc)(pmModel *model); +pmModelFunc pmModelFunc_GetFunction (pmModelType type); +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); +int pmModelParameterCount(pmModelType type); static void modelFree(pmModel *tmp) Index: /branches/rel10_ifa/psModules/src/objects/pmModel.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModel.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmModel.h (revision 6556) @@ -1,15 +1,11 @@ /** @file pmObjects.h * - * The process of finding, measuring, and classifying astronomical sources on - * images is one of the critical tasks of the IPP or any astronomical software - * system. This file will define structures and functions related to the task - * of source detection and measurement. The elements defined in this section - * are generally low-level components which can be connected together to - * construct a complete object measurement suite. + * Functions to define and manipulate object models * * @author GLG, MHPCC + * @author EAM, IfA * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -61,5 +57,5 @@ pmModel *pmModelCopy (pmModel *model); -psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row); /** pmModelAdd() Index: /branches/rel10_ifa/psModules/src/objects/pmModelGroup.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModelGroup.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmModelGroup.c (revision 6556) @@ -1,8 +1,35 @@ -# include "pmModelGroup.h" - +/** @file pmModelGroup.c + * + * Functions to define and manipulate object model attributes + * + * @author GLG, MHPCC + * @author EAM, IfA + * + * @version $Revision: 1.2.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" + +// XXX shouldn't these be defined for us in pslib.h ??? double hypot(double x, double y); double sqrt (double x); -#include "psEllipse.h" #include "models/pmModel_GAUSS.c" #include "models/pmModel_PGAUSS.c" @@ -159,2 +186,28 @@ return (models[type].name); } + +/****************************************************************************** + pmSourceModelGuess(source, model): This function allocates a new + pmModel structure based on the given modelType specified in the argument list. + The corresponding pmModelGuess function is returned, and used to + supply the values of the params array in the pmModel structure. + + XXX: Many parameters are based on the src->moments structure, which is in + image, not subImage coords. Therefore, the calls to the model evaluation + functions will be in image, not subImage coords. Remember this. +*****************************************************************************/ +pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); + modelGuessFunc(model, source); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(model); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmModelGroup.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModelGroup.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmModelGroup.h (revision 6556) @@ -9,16 +9,84 @@ * @author EAM, IfA * - * @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.2.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * */ -#include "pmObjects.h" -#include "pmPSF.h" -/** - * - * This function returns the number of parameters used by the listed function. - * + +# ifndef PM_MODEL_GROUP_H +# define PM_MODEL_GROUP_H + +// This function is the model chi-square minimization function for this model. +typedef psMinimizeLMChi2Func pmModelFunc; + +// This function returns the integrated flux for the given model parameters. +typedef psF64 (*pmModelFlux)(const psVector *params); + + +// This function returns the radius at which the given model and parameters +// achieves the given flux. +typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + +/* This function sets the model parameter limits vectors for the given model + */ +typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); + +/* This function provides the model guess parameters based on the details of + * the given source. + */ +typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + + +/* This function constructs the PSF model for the given source based on the + * supplied psf and the EXT model for the object. + */ +typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf); + +/* This function returns the success / failure status of the given model fit + */ +typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +/* Every model instance belongs to a class of models, defined by the value of + * the pmModelType type entry. Various functions need access to information about + * each of the models. Some of this information varies from model to model, and + * may depend on the current parameter values or other data quantities. In order + * to keep the code from requiring the information about each model to be coded + * into the low-level fitting routines, we define a collection of functions which + * allow us to abstract this type of model-dependent information. These generic + * functions take the model type and return the corresponding function pointer + * for the specified model. Each model is defined by creating this collection of + * specific functions, and placing them in a single file for each model. We + * define the following structure to carry the collection of information about + * the models. + */ +typedef struct +{ + char *name; + int nParams; + pmModelFunc modelFunc; + pmModelFlux modelFlux; + pmModelRadius modelRadius; + pmModelLimits modelLimits; + pmModelGuessFunc modelGuessFunc; + pmModelFromPSFFunc modelFromPSFFunc; + pmModelFitStatusFunc modelFitStatusFunc; +} +pmModelGroup; + +// allocate a pmModelGroup to hold nModels entries +pmModelGroup *pmModelGroupAlloc (int nModels); + +// initialize the internal (static) model group with the default models +void pmModelGroupInit (void); + +// free the internal (static) model group +void pmModelGroupCleanup (void); + +// add a new model to the internal (static) model group +void pmModelGroupAdd (pmModelGroup *model); + +/* This function returns the number of parameters used by the listed function. */ int pmModelParameterCount( @@ -27,8 +95,5 @@ -/** - * - * This function returns the user-space model names for the specified model type. - * +/* This function returns the user-space model names for the specified model type. */ char *pmModelGetType( @@ -46,62 +111,4 @@ ); - -#ifndef PM_MODEL_GROUP_H -#define PM_MODEL_GROUP_H - -/** - * - * This function is the model chi-square minimization function for this model. - * - */ -typedef psMinimizeLMChi2Func pmModelFunc; - - -/** - * - * This function returns the integrated flux for the given model parameters. - */ -typedef psF64 (*pmModelFlux)(const psVector *params); - - -/** - * - * This function returns the radius at which the given model and parameters - * achieves the given flux. - * - */ -typedef psF64 (*pmModelRadius)(const psVector *params, double flux); - -/** - * - * This function sets the model parameter limits vectors for the given model - * - */ -typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); - -/** - * - * This function provides the model guess parameters based on the details of - * the given source. - * - */ -typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); - - -/** - * - * This function constructs the PSF model for the given source based on the - * supplied psf and the EXT model for the object. - * - */ -typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf); - -/** - * - * This function returns the success / failure status of the given model fit - * - */ -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - /** * @@ -110,5 +117,4 @@ * */ - /** @@ -177,47 +183,17 @@ +/** pmSourceModelGuess() + * + * Convert available data to an initial guess for the given model. This + * function allocates a pmModel entry for the pmSource structure based on the + * provided model selection. The method of defining the model parameter guesses + * are specified for each model below. The guess values are placed in the model + * parameters. The function returns TRUE on success or FALSE on failure. + * + */ +pmModel *pmSourceModelGuess( + pmSource *source, ///< The input pmSource + pmModelType model ///< The type of model to be created. +); - -/** - * - * Every model instance belongs to a class of models, defined by the value of - * the pmModelType type entry. Various functions need access to information about - * each of the models. Some of this information varies from model to model, and - * may depend on the current parameter values or other data quantities. In order - * to keep the code from requiring the information about each model to be coded - * into the low-level fitting routines, we define a collection of functions which - * allow us to abstract this type of model-dependent information. These generic - * functions take the model type and return the corresponding function pointer - * for the specified model. Each model is defined by creating this collection of - * specific functions, and placing them in a single file for each model. We - * define the following structure to carry the collection of information about - * the models. - * - */ -typedef struct -{ - char *name; - int nParams; - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuessFunc; - pmModelFromPSFFunc modelFromPSFFunc; - pmModelFitStatusFunc modelFitStatusFunc; -} -pmModelGroup; - -// allocate a pmModelGroup to hold nModels entries -pmModelGroup *pmModelGroupAlloc (int nModels); - -// initialize the internal (static) model group with the default models -void pmModelGroupInit (void); - -// free the internal (static) model group -void pmModelGroupCleanup (void); - -// add a new model to the internal (static) model group -void pmModelGroupAdd (pmModelGroup *model); - -# endif +# endif /* PM_MODEL_GROUP_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmMoments.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmMoments.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmMoments.c (revision 6556) @@ -1,11 +1,11 @@ -/** @file pmObjects.c +/** @file pmMoments.c * - * This file will ... + * Functions defining the pmMoments structure * * @author GLG, MHPCC - * @author EAM, IfA: significant modifications. + * @author EAM, IfA * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,10 +13,6 @@ */ -#include -#include -#include #include "pslib.h" -#include "pmObjects.h" -#include "pmModelGroup.h" +#include "pmMoments.h" /****************************************************************************** @@ -42,177 +38,2 @@ } -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -pmSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - pmSource - pmSource->peak - pmSource->pixels - pmSource->weight - pmSource->mask - -XXX: The peak calculations are done in image coords, not subImage coords. - -XXX EAM : this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons -*****************************************************************************/ -# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -bool pmSourceMoments(pmSource *source, - psF32 radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - psF32 xPeak = source->peak->x; - psF32 yPeak = source->peak->y; - psF32 xOff = source->pixels->col0 - source->peak->x; - psF32 yOff = source->pixels->row0 - source->peak->y; - - // XXX why do I get different results for these two methods of finding Sx? - // XXX Sx, Sy would be better measured if we clip pixels close to sky - // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? - // We loop through all pixels in this subimage (source->pixels), and for each - // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; - psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if ((vMsk != NULL) && *vMsk) { - vMsk++; - continue; - } - - psF32 xDiff = col + xOff; - psF32 yDiff = row + yOff; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) { - if (vMsk != NULL) - vMsk++; - continue; - } - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : check for valid S/N in pixel - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) { - if (vMsk != NULL) - vMsk++; - continue; - } - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - XY += xDiff * yWght; - X2 += xDiff * xWght; - Y2 += yDiff * yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - if (vMsk != NULL) - vMsk++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return (false); - } - - psTrace (".psModules.pmSourceMoments", 5, - "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace (".psModules.pmSourceMoments", 3, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return (false); - } - - source->moments->x = x + xPeak; - source->moments->y = y + yPeak; - - // XXX EAM : Sxy needs to have x*y subtracted - source->moments->Sxy = XY/Sum - x*y; - source->moments->Sum = Sum; - source->moments->SN = Sum / sqrt(Var); - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - - // XXX EAM : these values can be negative, so we need to limit the range - source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); - source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); - - psTrace (".psModules.pmSourceMoments", 4, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, source->moments->Sxy); - - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - Index: /branches/rel10_ifa/psModules/src/objects/pmMoments.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmMoments.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmMoments.h (revision 6556) @@ -1,19 +1,17 @@ -/** @file pmObjects.h +/** @file pmMoments.h * - * The process of finding, measuring, and classifying astronomical sources on - * images is one of the critical tasks of the IPP or any astronomical software - * system. This file will define structures and functions related to the task - * of source detection and measurement. The elements defined in this section - * are generally low-level components which can be connected together to - * construct a complete object measurement suite. + * Definitions of the moments structure * * @author GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * */ + +# ifndef PM_MOMENTS_H +# define PM_MOMENTS_H /** pmMoments data structure @@ -45,16 +43,3 @@ pmMoments *pmMomentsAlloc(); -/** pmSourceMoments() - * - * Measure source moments for the given source, using the value of - * source.moments.sky provided as the local background value and the peak - * coordinates as the initial source location. The resulting moment values are - * applied to the source.moments entry, and the source is returned. The moments - * are measured within the given circular radius of the source.peak coordinates. - * The return value indicates the success (TRUE) of the operation. - * - */ -bool pmSourceMoments( - pmSource *source, ///< The input pmSource for which moments will be computed - float radius ///< Use a circle of pixels around the peak -); +# endif Index: /branches/rel10_ifa/psModules/src/objects/pmPSF.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSF.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmPSF.c (revision 6556) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.4.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,5 +18,11 @@ #include -#include "pmObjects.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmGrowthCurve.h" #include "pmPSF.h" #include "pmModelGroup.h" Index: /branches/rel10_ifa/psModules/src/objects/pmPSF.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSF.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmPSF.h (revision 6556) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1.22.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.1.22.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -15,6 +15,4 @@ # ifndef PM_PSF_H # define PM_PSF_H - -#include "pmGrowthCurve.h" /** pmPSF data structure Index: /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c (revision 6556) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.4.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,8 +13,16 @@ # include -# include "pmObjects.h" -# include "pmPSF.h" -# include "pmPSFtry.h" -# include "pmModelGroup.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmGrowthCurve.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmModelGroup.h" +#include "pmSourceFitModel.h" +#include "pmSourcePhotometry.h" // ******** pmPSFtry functions ************************************************** @@ -104,7 +112,7 @@ // fit model as EXT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED); // exclude the poor fits @@ -138,5 +146,5 @@ y = source->peak->y; - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); status = pmSourceFitModel (source, modelPSF, true); @@ -154,5 +162,9 @@ // XXX : use a different estimator for the local sky? // XXX : pass 'source' as input? - if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { + if (!pmSourcePhotometryModel (&fitMag, modelPSF)) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; + goto next_source; + } + if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) { psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; goto next_source; @@ -164,5 +176,5 @@ next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED); } Index: /branches/rel10_ifa/psModules/src/objects/pmPeaks.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPeaks.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmPeaks.c (revision 6556) @@ -1,11 +1,11 @@ -/** @file pmObjects.c +/** @file pmPeaks.c * - * This file will ... + * This file defines functions to detect and manipulate peaks in images * * @author GLG, MHPCC * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,6 +17,5 @@ #include #include "pslib.h" -#include "pmObjects.h" -#include "pmModelGroup.h" +#include "pmPeaks.h" /****************************************************************************** @@ -110,7 +109,6 @@ } -// XXX are these private or public? // psSort comparison function for peaks -int pmComparePeakAscend (const void **a, const void **b) +int pmPeaksCompareAscend (const void **a, const void **b) { psTrace(__func__, 3, "---- %s() begin ----\n", __func__); @@ -133,5 +131,5 @@ // psSort comparison function for peaks -int pmComparePeakDescend (const void **a, const void **b) +int pmPeaksCompareDescend (const void **a, const void **b) { psTrace(__func__, 3, "---- %s() begin ----\n", __func__); Index: /branches/rel10_ifa/psModules/src/objects/pmPeaks.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPeaks.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmPeaks.h (revision 6556) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -138,3 +138,6 @@ ); +int pmPeaksCompareAscend (const void **a, const void **b); +int pmPeaksCompareDescend (const void **a, const void **b); + # endif /* PM_PEAKS_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSource.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSource.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSource.c (revision 6556) @@ -1,11 +1,11 @@ -/** @file pmObjects.c +/** @file pmSource.c * - * This file will ... + * Functions to define and manipulate sources on images * * @author GLG, MHPCC * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,6 +17,10 @@ #include #include "pslib.h" -#include "pmObjects.h" -#include "pmModelGroup.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" static void sourceFree(pmSource *tmp) @@ -50,6 +54,6 @@ tmp->modelPSF = NULL; tmp->modelEXT = NULL; - tmp->type = PM_SOURCE_UNKNOWN; - tmp->mode = PM_SOURCE_DEFAULT; + tmp->type = PM_SOURCE_TYPE_UNKNOWN; + tmp->mode = PM_SOURCE_MODE_DEFAULT; psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); @@ -219,5 +223,5 @@ // XXX EAM : this lets us takes the single highest peak - psArraySort (peaks, pmComparePeakDescend); + psArraySort (peaks, pmPeaksCompareDescend); clump = peaks->data[0]; psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); @@ -331,5 +335,5 @@ // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2); - int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PM_SOURCE_MASK_SATURATED); // saturated star (size consistent with PSF or larger) @@ -338,6 +342,6 @@ big = true; if ((Nsatpix > 1) && big) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_SATSTAR; + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_SATSTAR; Nsatstar ++; continue; @@ -346,6 +350,6 @@ // saturated object (not a star, eg bleed trails, hot pixels) if (Nsatpix > 1) { - tmpSrc->type = PM_SOURCE_SATURATED; - tmpSrc->mode = PM_SOURCE_DEFAULT; + tmpSrc->type = PM_SOURCE_TYPE_SATURATED; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; Nsat ++; continue; @@ -356,6 +360,6 @@ // only set candidate defects if if ((sigX < 0.05) || (sigY < 0.05)) { - tmpSrc->type = PM_SOURCE_DEFECT; - tmpSrc->mode = PM_SOURCE_DEFAULT; + tmpSrc->type = PM_SOURCE_TYPE_DEFECT; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; Ncr ++; continue; @@ -364,6 +368,6 @@ // likely unsaturated extended source (too large to be stellar) if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - tmpSrc->type = PM_SOURCE_EXTENDED; - tmpSrc->mode = PM_SOURCE_DEFAULT; + tmpSrc->type = PM_SOURCE_TYPE_EXTENDED; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; Next ++; continue; @@ -378,6 +382,6 @@ psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_PSFSTAR; + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_PSFSTAR; Npsf ++; continue; @@ -385,6 +389,6 @@ // random type of star - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_DEFAULT; + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; } @@ -409,2 +413,191 @@ } +/****************************************************************************** +pmSourceMoments(source, radius): this function takes a subImage defined in the +pmSource data structure, along with the peak location, and determines the +various moments associated with that peak. + +Requires the following to have been created: + pmSource + pmSource->peak + pmSource->pixels + pmSource->weight + pmSource->mask + +XXX: The peak calculations are done in image coords, not subImage coords. + +XXX EAM : this version clips input pixels on S/N +XXX EAM : this version returns false for several reasons +*****************************************************************************/ +# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + +bool pmSourceMoments(pmSource *source, + psF32 radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); + + // + // XXX: Verify the setting for sky if source->moments == NULL. + // + psF32 sky = 0.0; + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } else { + sky = source->moments->Sky; + } + + // + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // X2 = SUM (x - xc)^2 * (z - sky) + // XY = SUM (x - xc)*(y - yc)*(z - sky) + // + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + psF32 Sum = 0.0; + psF32 Var = 0.0; + psF32 X1 = 0.0; + psF32 Y1 = 0.0; + psF32 X2 = 0.0; + psF32 Y2 = 0.0; + psF32 XY = 0.0; + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + psF32 xOff = source->pixels->col0 - source->peak->x; + psF32 yOff = source->pixels->row0 - source->peak->y; + + // XXX why do I get different results for these two methods of finding Sx? + // XXX Sx, Sy would be better measured if we clip pixels close to sky + // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? + // We loop through all pixels in this subimage (source->pixels), and for each + // pixel that is not masked, AND within the radius of the peak pixel, we + // proceed with the moments calculation. need to do two loops for a + // numerically stable result. first loop: get the sums. + // XXX EAM : mask == 0 is valid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->weight->data.F32[row]; + psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if ((vMsk != NULL) && *vMsk) { + vMsk++; + continue; + } + + psF32 xDiff = col + xOff; + psF32 yDiff = row + yOff; + + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) { + if (vMsk != NULL) + vMsk++; + continue; + } + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (PS_SQR(pDiff) < wDiff) { + if (vMsk != NULL) + vMsk++; + continue; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + XY += xDiff * yWght; + X2 += xDiff * xWght; + Y2 += yDiff * yWght; + + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + if (vMsk != NULL) + vMsk++; + } + } + + // if we have less than (1/4) of the possible pixels, force a retry + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 0.75*R2) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); + + // + // first moment X = X1/Sum + xc + // second moment X = sqrt (X2/Sum - (X1/Sum)^2) + // Sxy = XY / Sum + // + x = X1/Sum; + y = Y1/Sum; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 3, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + // XXX EAM : Sxy needs to have x*y subtracted + source->moments->Sxy = XY/Sum - x*y; + source->moments->Sum = Sum; + source->moments->SN = Sum / sqrt(Var); + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + // XXX EAM : these values can be negative, so we need to limit the range + source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); + source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + +pmModel *pmSourceSelectModel (pmSource *source) +{ + switch (source->type) { + case PM_SOURCE_TYPE_STAR: + return source->modelPSF; + + case PM_SOURCE_TYPE_EXTENDED: + return source->modelEXT; + + default: + return NULL; + } + return NULL; +} Index: /branches/rel10_ifa/psModules/src/objects/pmSource.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSource.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSource.h (revision 6556) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -81,5 +81,7 @@ pmSourceMode mode; ///< Best identification of object. psArray *blends; - float fitMag; + float psfMag; + float extMag; + float errMag; float apMag; psRegion region; // area on image covered by selected pixels @@ -125,15 +127,19 @@ * */ -// XXX: Uncommenting the pmReadout causes compile errors. bool pmSourceDefinePixels( - pmSource *mySource, ///< Add comment. - pmReadout *readout, ///< Add comment. - psF32 x, ///< Add comment. - psF32 y, ///< Add comment. - psF32 Radius ///< Add comment. -); - -bool pmSourceDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); -bool pmSourceRedefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); + pmSource *mySource, ///< source to be re-defined + const pmReadout *readout, ///< base the source on this readout + psF32 x, ///< center coords of source + psF32 y, ///< center coords of source + psF32 Radius ///< size of box on source +); + +bool pmSourceRedefinePixels ( + pmSource *mySource, ///< source to be re-defined + const pmReadout *readout, ///< base the source on this readout + psF32 x, ///< center coords of source + psF32 y, ///< center coords of source + psF32 Radius ///< size of box on source +); /** pmSourcePSFClump() @@ -188,3 +194,22 @@ +/** pmSourceMoments() + * + * Measure source moments for the given source, using the value of + * source.moments.sky provided as the local background value and the peak + * coordinates as the initial source location. The resulting moment values are + * applied to the source.moments entry, and the source is returned. The moments + * are measured within the given circular radius of the source.peak coordinates. + * The return value indicates the success (TRUE) of the operation. + * + */ +bool pmSourceMoments( + pmSource *source, ///< The input pmSource for which moments will be computed + float radius ///< Use a circle of pixels around the peak +); + + +// select the model used for this source +pmModel *pmSourceSelectModel (pmSource *source); + # endif /* PM_SOURCE_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c (revision 6556) @@ -1,3 +1,27 @@ -# include "psphot.h" +/** @file pmSourceContour.c + * + * Functions to measure the local sky and sky variance for sources on images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceContour.h" /****************************************************************************** Index: /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h (revision 6556) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,25 +13,25 @@ # define PM_SOURCE_CONTOUR_H -psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold) -{ +psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold); - /** pmSourceContour() - * - * Find points in a contour for the given source at the given level. If type - * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, - * running along each pixel row until the level is crossed, then interpolating to - * the level coordinate for that row. This is done for each row, with the - * starting point determined by the midpoint of the previous row, until the - * starting point has a value below the contour level. The returned contour - * consists of two vectors giving the x and y coordinates of the contour levels. - * This function may be used as part of the model guess inputs. Other contour - * types may be specified in the future for more refined contours (TBD) - * - */ - psArray *pmSourceContour_Crude( - pmSource *source, ///< The input pmSource - const psImage *image, ///< The input image (float) (this arg should be removed) - float level ///< The level of the contour - ); - # endif /* PM_SOURCE_PHOTOMETRY_H */ +/** pmSourceContour() + * + * Find points in a contour for the given source at the given level. If type + * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, + * running along each pixel row until the level is crossed, then interpolating to + * the level coordinate for that row. This is done for each row, with the + * starting point determined by the midpoint of the previous row, until the + * starting point has a value below the contour level. The returned contour + * consists of two vectors giving the x and y coordinates of the contour levels. + * This function may be used as part of the model guess inputs. Other contour + * types may be specified in the future for more refined contours (TBD) + * + */ +psArray *pmSourceContour_Crude( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) (this arg should be removed) + float level ///< The level of the contour +); + +# endif /* PM_SOURCE_PHOTOMETRY_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6556) @@ -1,11 +1,11 @@ -/** @file pmObjects.c +/** @file pmSourceFitModel.c * - * This file will ... + * fit single source models to image pixels * + * @author EAM, IfA * @author GLG, MHPCC - * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,6 +17,14 @@ #include #include "pslib.h" -#include "pmObjects.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" #include "pmModelGroup.h" +#include "pmSourceFitModel.h" // save a static values so they may be set externally @@ -169,5 +177,5 @@ } - source->mode |= PM_SOURCE_FITTED; + source->mode |= PM_SOURCE_MODE_FITTED; psFree(x); @@ -186,43 +194,2 @@ return(rc); } - -pmModel *pmSourceSelectModel (pmSource *source) -{ - switch (source->type) { - case PM_SOURCE_STAR: - return source->modelPSF; - - case PM_SOURCE_EXTENDED: - return source->modelEXT; - - default: - return NULL; - } - return NULL; -} - -/****************************************************************************** - pmSourceModelGuess(source, model): This function allocates a new - pmModel structure based on the given modelType specified in the argument list. - The corresponding pmModelGuess function is returned, and used to - supply the values of the params array in the pmModel structure. - - XXX: Many parameters are based on the src->moments structure, which is in - image, not subImage coords. Therefore, the calls to the model evaluation - functions will be in image, not subImage coords. Remember this. -*****************************************************************************/ -pmModel *pmSourceModelGuess(pmSource *source, - pmModelType modelType) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - - pmModel *model = pmModelAlloc(modelType); - - pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); - modelGuessFunc(model, source); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(model); -} - Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6556) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -33,20 +33,3 @@ ); - -pmModel *pmSourceSelectModel (pmSource *source); - -/** pmSourceModelGuess() - * - * Convert available data to an initial guess for the given model. This - * function allocates a pmModel entry for the pmSource structure based on the - * provided model selection. The method of defining the model parameter guesses - * are specified for each model below. The guess values are placed in the model - * parameters. The function returns TRUE on success or FALSE on failure. - * - */ -pmModel *pmSourceModelGuess( - pmSource *source, ///< The input pmSource - pmModelType model ///< The type of model to be created. -); - # endif /* PM_SOURCE_FIT_MODEL_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6556) @@ -1,8 +1,8 @@ -/** @file pmSourceFitSet.h +/** @file pmSourceFitSet.c * * @author EAM, IfA * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -10,11 +10,21 @@ */ +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceFitSet.h" + // sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, // nPar = nSrc*(nOnePar - 1) + 1 - -# include -# include -# include -# include "pmSource.h" static pmModelFunc mFunc; @@ -110,5 +120,5 @@ // set the static variables - pmSourceFitSetInit (model->type); + pmModelFitSetInit (model->type); int nSrc = modelSet->n; @@ -272,5 +282,5 @@ } - source->mode |= PM_SOURCE_FITTED; + source->mode |= PM_SOURCE_MODE_FITTED; psFree(x); Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c (revision 6556) @@ -1,2 +1,27 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" // translations between psphot object types and dophot object types @@ -6,18 +31,18 @@ switch (source->type) { - case PM_SOURCE_DEFECT: - case PM_SOURCE_SATURATED: + case PM_SOURCE_TYPE_DEFECT: + case PM_SOURCE_TYPE_SATURATED: return (8); - case PM_SOURCE_STAR: - if (source->mode & PM_SOURCE_SATSTAR) + case PM_SOURCE_TYPE_STAR: + if (source->mode & PM_SOURCE_MODE_SATSTAR) return (10); - if (source->mode & PM_SOURCE_POOR) + if (source->mode & PM_SOURCE_MODE_POOR) return (7); - if (source->mode & PM_SOURCE_FAIL) + if (source->mode & PM_SOURCE_MODE_FAIL) return (4); return (1); - case PM_SOURCE_EXTENDED: + case PM_SOURCE_TYPE_EXTENDED: return (2); @@ -29,6 +54,6 @@ /***** Text Output Methods *****/ - -bool pmSourcesWriteText (psArray *sources, char *filename) +# if (0) + bool pmSourcesWriteText (psArray *sources, char *filename) { @@ -50,3 +75,3 @@ return true; } - +# endif Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO.h (revision 6556) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO.h (revision 6556) @@ -0,0 +1,26 @@ +/** @file pmSourceIO.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_IO_H +# define PM_SOURCE_IO_H + +int pmSourceDophotType (pmSource *source); +bool pmSourcesWriteText (psArray *sources, char *filename); + +bool pmSourcesWriteOBJ (psArray *sources, char *filename); +bool pmSourcesWriteSX (psArray *sources, char *filename); +bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header); +bool pmSourcesWriteCMF (psFits *fits, psArray *sources, psMetadata *header, char *extname); + +psArray *pmSourcesReadCMP (char *filename, psMetadata *header); +psArray *pmSourcesReadCMF (psFits *fits, psMetadata *header); + +# endif /* PM_SOURCE_IO_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c (revision 6556) @@ -1,3 +1,28 @@ -# include "pmSource.h" +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" // elixir-style FITS table output (header + table in 1st extension) @@ -12,6 +37,8 @@ int i, type; psF32 *PAR, *dPAR; - float dmag, lsky; + float lsky; bool status; + EllipseShape shape; + EllipseAxes axes; // find config information for output header @@ -23,5 +50,5 @@ for (i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); + pmModel *model = pmSourceSelectModel (source); if (model == NULL) continue; @@ -30,20 +57,24 @@ dPAR = model->dparams->data.F32; - dmag = dPAR[1] / PAR[1]; type = pmSourceDophotType (source); lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = EllipseShapeToAxes (shape); row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", 32.0); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", source->apMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", PS_MIN (99.0, source->psfMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", PS_MIN (99.0, source->extMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", PS_MIN (99.0, source->apMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->errMag)); psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky); - psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", type); - psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); - psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); - psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0"); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_U8, "", type); psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123"); @@ -62,11 +93,25 @@ // read in a readout from the fits file -psArray *pmReadoutReadCMF (psFits *fits, psMetadata *header) +psArray *pmSourcesReadCMF (psFits *fits, psMetadata *header) { bool status; + psF32 *PAR, *dPAR; + EllipseShape shape; + EllipseAxes axes; + float lsky; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } // find config information for output header float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; psArray *table = psFitsReadTable (fits); @@ -76,19 +121,33 @@ // convert the table to the pmSource entries + // XXX need to chooose PSF vs EXT, based on type? for (int i = 0; i < table->n; i++) { - pmStar *source = pmSourceAlloc (); + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; psMetadata *row = table->data[i]; - source->x = psMetadataLookupF32 (&status, row, "X_PIX"); - source->y = psMetadataLookupF32 (&status, row, "Y_PIX"); - source->fitMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; - source->galMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; + lsky = psMetadataLookupF32 (&status, row, "LOG_SKY"); + PAR[0] = pow(lsky, 10.0); + PAR[2] = psMetadataLookupF32 (&status, row, "X_PIX"); + PAR[3] = psMetadataLookupF32 (&status, row, "Y_PIX"); + axes.major = psMetadataLookupF32 (&status, row, "FWHM_X"); + axes.minor = psMetadataLookupF32 (&status, row, "FWHM_Y"); + axes.theta = psMetadataLookupF32 (&status, row, "THETA"); + + shape = EllipseAxesToShape (axes); + + PAR[4] = shape.sx; + PAR[5] = shape.sy; + PAR[6] = shape.sxy; + + source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; + source->extMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; + source->errMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; - source->dMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); - source->sky = psMetadataLookupF32 (&status, row, "LOG_SKY"); - source->sx = psMetadataLookupF32 (&status, row, "FWHM_X"); - source->sy = psMetadataLookupF32 (&status, row, "FWHM_Y"); - source->theta = psMetadataLookupF32 (&status, row, "THETA"); sources->data[i] = source; Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c (revision 6556) @@ -1,3 +1,29 @@ -# include "pmSource.h" +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" // elixir-style pseudo FITS table (header + ascii list) @@ -10,4 +36,6 @@ float dmag, lsky; bool status; + EllipseShape shape; + EllipseAxes axes; // find config information for output header @@ -19,5 +47,5 @@ for (i = nSrc = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); + pmModel *model = pmSourceSelectModel (source); if (model == NULL) continue; @@ -49,5 +77,5 @@ for (i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); + pmModel *model = pmSourceSelectModel (source); if (model == NULL) continue; @@ -56,20 +84,24 @@ dPAR = model->dparams->data.F32; - dmag = dPAR[1] / PAR[1]; type = pmSourceDophotType (source); lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = EllipseShapeToAxes (shape); + psLineInit (line); - psLineAdd (line, "%6.1f ", PAR[2]); - psLineAdd (line, "%6.1f ", PAR[3]); - psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); - psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); - psLineAdd (line, "%2d ", type); - psLineAdd (line, "%3.1f ", lsky); - psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal - psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); - psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' - psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' - psLineAdd (line, "%5.1f\n", 0); // should be theta + psLineAdd (line, "%6.1f ", PAR[2]); + psLineAdd (line, "%6.1f ", PAR[3]); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->psfMag + ZERO_POINT)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->extMag + ZERO_POINT)); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->apMag + ZERO_POINT)); + psLineAdd (line, "%6.2f ", axes.major); + psLineAdd (line, "%6.2f ", axes.minor); + psLineAdd (line, "%5.1f\n", axes.theta); fwrite (line->line, 1, line->Nline, f); } @@ -82,12 +114,20 @@ // elixir-style pseudo FITS table (header + ascii list) -psArray *sources pmSourcesReadCMP (char *filename, psMetadata *header) +psArray *pmSourcesReadCMP (char *filename, psMetadata *header) { - int i, type, nSrc; - psMetadataItem *mdi; + int Ninstar; psF32 *PAR, *dPAR; - float dmag, lsky; bool status; + EllipseShape shape; + EllipseAxes axes; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } // find config information for output header @@ -113,10 +153,10 @@ // prepare array to store data - int nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS"); + int nStars = psMetadataLookupS32 (&status, header, "NSTARS"); psArray *sources = psArrayAlloc (nStars); sources->n = 0; // we have fixed bytes / line : use that info - char *buffer = psAlloc (NBYTES_STAR*nStars); + char *buffer = psAlloc (BYTES_STAR*nStars); int Nextra = 0; @@ -161,24 +201,36 @@ /* extract data for stars */ - int Ninstar = nbytes / BYTES_STAR; + Ninstar = nbytes / BYTES_STAR; Nextra = nbytes % BYTES_STAR; - for (j = 0; j < Ninstar; j++, N++) { - psString line = psStringNCopy (&buffer[j*BYTES_STAR]); + for (int j = 0; j < Ninstar; j++) { + psString line = psStringNCopy (&buffer[j*BYTES_STAR], BYTES_STAR); psList *list = psStringSplit (line, " "); psArray *array = psListToArray (list); pmSource *source = pmSourceAlloc (); - - source->x = atof (array->data[0]); - source->y = atof (array->data[1]); - source->fitMag = atof (array->data[2]); - source->dMag = 0.001*atof (array->data[3]); - source->type = atof (array->data[4]); - source->sky = pow (atof (array->data[5]), 10.0); - source->galMag = atof (array->data[6]); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + PAR[0] = pow (atof (array->data[5]), 10.0); + PAR[2] = atof (array->data[0]); + PAR[3] = atof (array->data[1]); + source->psfMag = atof (array->data[2]); + source->extMag = atof (array->data[6]); + source->errMag = atof (array->data[3]) / 1000.0; source->apMag = atof (array->data[7]); - source->sx = atof (array->data[8]); - source->sy = atof (array->data[9]); - source->theta = atof (array->data[10]); + axes.major = atof (array->data[8]); + axes.minor = atof (array->data[9]); + axes.theta = atof (array->data[10]); + + shape = EllipseAxesToShape (axes); + + PAR[4] = shape.sx; + PAR[5] = shape.sy; + PAR[6] = shape.sxy; + + // source->type = atof (array->data[4]); psArrayAdd (sources, 100, source); Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c (revision 6556) @@ -1,3 +1,29 @@ -# include "pmSource.h" +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" // dophot-style output list with fixed line width @@ -8,4 +34,6 @@ psF32 *PAR, *dPAR; float dmag, apResid; + EllipseShape shape; + EllipseAxes axes; psTimerStart ("string"); @@ -22,5 +50,5 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); + pmModel *model = pmSourceSelectModel (source); if (model == NULL) continue; @@ -31,5 +59,14 @@ dmag = dPAR[1] / PAR[1]; type = pmSourceDophotType (source); - apResid = source->apMag - source->fitMag; + if ((source->apMag < 99.0) && (source->psfMag < 99.0)) { + apResid = source->apMag - source->psfMag; + } else { + apResid = 0.0; + } + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = EllipseShapeToAxes (shape); psLineInit (line); @@ -37,11 +74,11 @@ psLineAdd (line, "%8.2f", PAR[2]); psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", source->fitMag); + psLineAdd (line, "%8.3f", source->psfMag); psLineAdd (line, "%6.3f", dmag); psLineAdd (line, "%9.2f", PAR[0]); - psLineAdd (line, "%9.3f", PAR[4]); - psLineAdd (line, "%9.3f", PAR[5]); - psLineAdd (line, "%7.2f", PAR[6]); - psLineAdd (line, "%8.3f", 99.999); + psLineAdd (line, "%9.3f", axes.major); + psLineAdd (line, "%9.3f", axes.minor); + psLineAdd (line, "%7.2f", axes.theta); + psLineAdd (line, "%8.3f", source->extMag); psLineAdd (line, "%8.3f", source->apMag); psLineAdd (line, "%8.2f\n", apResid); @@ -50,50 +87,5 @@ fclose (f); psFree (line); - fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n); - - psTimerStart ("string"); - - f = fopen ("test.obj", "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj"); - return false; - } - - char *string; - // write sources with models - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) - continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - dmag = dPAR[1] / PAR[1]; - type = pmSourceDophotType (source); - apResid = source->apMag - source->fitMag; - - string = NULL; - psStringAppend (&string, "%3d", type); - psStringAppend (&string, "%8.2f", PAR[2]); - psStringAppend (&string, "%8.2f", PAR[3]); - psStringAppend (&string, "%8.3f", source->fitMag); - psStringAppend (&string, "%6.3f", dmag); - psStringAppend (&string, "%9.2f", PAR[0]); - psStringAppend (&string, "%9.3f", PAR[4]); - psStringAppend (&string, "%9.3f", PAR[5]); - psStringAppend (&string, "%7.2f", PAR[6]); - psStringAppend (&string, "%8.3f", 99.999); - psStringAppend (&string, "%8.3f", source->apMag); - psStringAppend (&string, "%8.2f\n", apResid); - fwrite (string, 1, strlen(string), f); - psFree (string); - } - fclose (f); - fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n); - + fprintf (stderr, "%f seconds for %d objects\n", psTimerMark ("string"), (int)sources->n); return true; } - -// XXX should we use psStringAppend or psLineAdd? Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c (revision 6556) @@ -1,3 +1,29 @@ -# include "pmSource.h" +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" // elixir-mode / sextractor-style output list with fixed line width @@ -6,5 +32,6 @@ psF32 *PAR, *dPAR; - float dmag; + EllipseShape shape; + EllipseAxes axes; psLine *line = psLineAlloc (110); // 110 is sextractor line length @@ -19,5 +46,5 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); + pmModel *model = pmSourceSelectModel (source); if (model == NULL) continue; @@ -28,19 +55,22 @@ // pmSourceSextractType (source, &type, &flags); - dmag = dPAR[1] / PAR[1]; + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = EllipseShapeToAxes (shape); psLineInit (line); - psLineAdd (line, "%5.2f", 0.0); // should be type + psLineAdd (line, "%5.2f", 0.0); // should be type psLineAdd (line, "%11.3f", PAR[2]); psLineAdd (line, "%11.3f", PAR[3]); - psLineAdd (line, "%9.4f", source->fitMag); - psLineAdd (line, "%9.4f", dmag); + psLineAdd (line, "%9.4f", source->psfMag); + psLineAdd (line, "%9.4f", source->errMag); psLineAdd (line, "%13.4f", PAR[0]); - psLineAdd (line, "%9.2f", 0.0); // should be FWHMx - psLineAdd (line, "%9.2f", 0.0); // should be FWHMy - psLineAdd (line, "%6.1f", 0.0); // should be Theta - psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO - psLineAdd (line, "%9.4f", source->apMag); - psLineAdd (line, "%4d\n", 0); // should be flags + psLineAdd (line, "%9.2f", axes.major); + psLineAdd (line, "%9.2f", axes.minor); + psLineAdd (line, "%6.1f", axes.theta); + psLineAdd (line, "%9.4f", source->extMag); + psLineAdd (line, "%9.4f", source->apMag); + psLineAdd (line, "%4d\n", 0); // should be flags fwrite (line->line, 1, line->Nline, f); } Index: /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (revision 6556) @@ -1,8 +1,8 @@ -/** @file pmSourcePhotometry.C +/** @file pmSourcePhotometry.c * * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -10,7 +10,49 @@ */ -// XXX EAM : the apMag should only be calculated for the brighter sources? -// XXX EAM : SN limit set by user? +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourcePhotometry.h" + +static float AP_MIN_SN = 0.0; + +bool pmSourceMagnitudesInit (psMetadata *config) +{ + + bool status; + + float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN"); + if (status) { + AP_MIN_SN = limit; + } + return true; +} + +/** + this function is used to calculate the three defined source magnitudes: + - apMag : only if S/N > AP_MIN_SN + : is optionally corrected for curve-of-growth if: + - the source is a STAR (PSF) + - the option is selected (how??) + - psfMag : all sources with non-NULL modelPSF + : is optionally corrected for aperture residual if: + - the source is a STAR (PSF) + - the option is selected (how??) + + - extMag : all sources with non-NULL modelEXT +**/ + // XXX EAM : masked region should be (optionally) elliptical +// XXX curve of growth is corrected to pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) { @@ -21,8 +63,9 @@ float rflux; float radius; + float SN; pmModel *model; switch (source->type) { - case PM_SOURCE_STAR: + case PM_SOURCE_TYPE_STAR: model = source->modelPSF; if (model == NULL) @@ -32,5 +75,5 @@ break; - case PM_SOURCE_EXTENDED: + case PM_SOURCE_TYPE_EXTENDED: model = source->modelEXT; if (model == NULL) @@ -44,29 +87,46 @@ } + if (model->dparams->data.F32[1] > 0) { + SN = model->params->data.F32[1] / model->dparams->data.F32[1]; + source->errMag = 1.0 / SN; + } else { + SN = 0.0; + source->errMag = 0.0; + } x = model->params->data.F32[2]; y = model->params->data.F32[3]; // replace source flux + // XXX test to see if source has been subtracted? pmModelAdd (source->pixels, source->mask, model, false, false); - // set aperture mask circle of PSF_FIT_RADIUS - psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED); - - // measure object photometry - status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask); - - // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + // set aperture mask circle to scaled radius + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PM_SOURCE_MASK_MARKED); + + // measure object model photometry + status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + + // measure object aperture photometry + if (SN > AP_MIN_SN) { + status = pmSourcePhotometryAper (&source->apMag, model, source->pixels, source->mask); + } else { + source->apMag = 99.0; + } + + // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS if (isPSF && (psf != NULL)) { - if (psf->growth != NULL) { - source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius); - } - - rflux = pow (10.0, 0.4*source->fitMag); - source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; - source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); + if (SN > AP_MIN_SN) { + if (psf->growth != NULL) { + source->apMag += pmGrowthCurveCorrect (psf->growth, radius); + } + rflux = pow (10.0, 0.4*source->psfMag); + source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; + } + source->psfMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); } // unmask aperture - psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PM_SOURCE_MASK_MARKED); // subtract object, leave local sky @@ -85,6 +145,58 @@ */ -// XXX split into ap, psf, ext mags? -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) +// return source model magnitude +bool pmSourcePhotometryModel (float *fitMag, pmModel *model) +{ + + float fitSum = 0; + *fitMag = 99.0; + + if (model == NULL) { + return false; + } + + // measure fitMag + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + if (fitSum <= 0) + return false; + if (!isfinite(fitSum)) + return false; + *fitMag = -2.5*log10(fitSum); + + return (true); +} + +// return source aperture magnitude +bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask) +{ + + float apSum = 0; + *apMag = 99.0; + + if (model == NULL) { + return false; + } + + float sky = model->params->data.F32[0]; + + // measure apMag + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) + continue; + apSum += image->data.F32[iy][ix] - sky; + } + } + if (apSum <= 0) + return false; + *apMag = -2.5*log10(apSum); + + return (true); +} + +# if (0) + // XXX split into ap, psf, ext mags? + bool pmSourcePhotometry (float *fitMag, float *apMag, pmModel *model, psImage *image, psImage *mask) { @@ -96,4 +208,5 @@ *obsMag = 99.0; + // measure fitMag pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); fitSum = modelFluxFunc (model->params); @@ -104,4 +217,5 @@ *fitMag = -2.5*log10(fitSum); + // measure apMag for (int ix = 0; ix < image->numCols; ix++) { for (int iy = 0; iy < image->numRows; iy++) { @@ -117,4 +231,5 @@ return (true); } +# endif float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) Index: /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h (revision 6556) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -25,7 +25,11 @@ * */ -bool pmSourcePhotometry( +bool pmSourcePhotometryModel( float *fitMag, ///< integrated fit magnitude - float *obsMag, ///< aperture flux magnitude + pmModel *model ///< model used for photometry +); + +bool pmSourcePhotometryAper( + float *apMag, ///< aperture flux magnitude pmModel *model, ///< model used for photometry psImage *image, ///< image pixels to be used Index: /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c (revision 6555) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c (revision 6556) @@ -1,40 +1,36 @@ +/** @file pmSourceSky.c + * + * Functions to measure the local sky and sky variance for sources on images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-09 03:14:23 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceSky.h" /****************************************************************************** -pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this -routine creates a new pmSource data structure and sets the following members: - ->pmPeak - ->pmMoments->sky +pmSource *pmSourceLocalSky(source, statsOptions, Radius): this +routine creates a new pmSource.moments element if needed and sets pmSource.pmMoments.sky The sky value is set from the pixels in the square annulus surrounding the peak pixel. -We simply create a subSet image and mask the inner pixels, then call -psImageStats on that subImage+mask. - -XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. - -XXX: Use static data structures for: - subImage - subImageMask - myStats - -XXX: ensure that the inner and out radius fit in the actual image. Should - we generate an error, or warning? Currently an error. - -XXX: Sync with IfA on whether the peak x/y coords are data structure coords, - or they use the image row/column offsets. -XXX EAM : peak->x,y uses parent coordinates - -XXX: Should we simply set pmSource->peak = peak? If so, should we increase -the reference counter? Or, should we copy the data structure? - -XXX: Currently the subimage always has an even number of rows/columns. Is - this correct? Since there is a center pixel, maybe it should have an - odd number of rows/columns. - -XXX: Use psTrace() for the print statements. - -XXX: Don't use separate structs for the subimage and mask. Use the source-> - members. +The source.pixels and source.mask must already exist *****************************************************************************/ @@ -60,8 +56,8 @@ srcRegion = psRegionForImage(mask, srcRegion); - psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); + psImageMaskRegion(mask, srcRegion, "OR", PM_SOURCE_MASK_MARKED); psStats *myStats = psStatsAlloc(statsOptions); myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); + psImageMaskRegion(mask, srcRegion, "AND", ~PM_SOURCE_MASK_MARKED); psF64 tmpF64; @@ -103,8 +99,8 @@ srcRegion = psRegionForImage(mask, srcRegion); - psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); + psImageMaskRegion(mask, srcRegion, "OR", PM_SOURCE_MASK_MARKED); psStats *myStats = psStatsAlloc(statsOptions); myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); + psImageMaskRegion(mask, srcRegion, "AND", ~PM_SOURCE_MASK_MARKED); psF64 tmpF64;