Changeset 6571
- Timestamp:
- Mar 11, 2006, 5:27:13 PM (20 years ago)
- Location:
- trunk/psphot/src
- Files:
-
- 2 added
- 5 deleted
- 18 edited
-
Makefile.am (modified) (3 diffs)
-
pmCellSetMask.c (deleted)
-
pmModelFitSet.c (deleted)
-
pmSourceContour.c (deleted)
-
pmSourceFitSet.c (deleted)
-
psModulesUtils.c (modified) (1 diff)
-
psModulesUtils.h (modified) (1 diff)
-
psphot.c (modified) (1 diff)
-
psphot.h (modified) (9 diffs)
-
psphotArguments.c (modified) (4 diffs)
-
psphotChoosePSF.c (modified) (1 diff)
-
psphotCleanup.c (modified) (1 diff)
-
psphotEnsemblePSF.c (modified) (1 diff)
-
psphotFitSet.c (added)
-
psphotImageMedian.c (modified) (6 diffs)
-
psphotMagnitudes.c (modified) (2 diffs)
-
psphotMaskCell.c (added)
-
psphotModelGroupInit.c (modified) (2 diffs)
-
psphotModelTest.c (modified) (7 diffs)
-
psphotOutput.c (modified) (3 diffs)
-
psphotParseCamera.c (modified) (2 diffs)
-
psphotReadout.c (modified) (5 diffs)
-
psphotSetup.c (deleted)
-
psphotSkyReplace.c (modified) (2 diffs)
-
testPSphotLoop.c (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/Makefile.am
r6481 r6571 12 12 13 13 psphot_SOURCES = \ 14 psphot.c \ 15 psphotModelGroupInit.c \ 16 psphotArguments.c \ 17 psphotParseCamera.c \ 18 testPSphotLoop.c \ 19 psphotMaskCell.c \ 14 20 psphotReadout.c \ 15 psphotImageStats.c \ 16 psphotImageBackground.c \ 21 psphotImageMedian.c \ 17 22 psphotFindPeaks.c \ 18 23 psphotSourceStats.c \ … … 24 29 psphotReplaceUnfit.c \ 25 30 psphotApResid.c \ 26 psphotOutput.c \ 27 psphotLoadPixels.c \ 28 psphotModelGroupInit.c \ 29 psphotGrowthCurve.c \ 30 psphotCleanup.c \ 31 psphotMagnitudes.c \ 32 psphotSkyReplace.c \ 31 33 psphotEvalPSF.c \ 32 34 psphotEvalFLT.c \ 33 psphotSourceFits.c \ 35 psphotSourceFits.c 36 37 junk = psphotGrowthCurve.c \ 38 psphotCleanup.c \ 34 39 psphotSortBySN.c \ 35 40 psphotDefinePixels.c \ 36 psphotMagnitudes.c \37 41 psphotRadiusChecks.c \ 38 psphotImageMedian.c \39 psphotSkyReplace.c \40 psLine.c \41 psSparse.c \42 psPolynomialUtils.c \43 42 psModulesUtils.c \ 44 43 pmSourceContour.c \ 45 44 pmModelFitSet.c \ 46 pmCellSetMask.c \47 45 pmSourceFitSet.c \ 48 psphot.c \ 49 psphotArguments.c \ 50 psphotParseCamera.c \ 51 psphotModelTest.c \ 52 psphotImageLoop.c 46 psphotModelTest.c 53 47 54 48 noinst_HEADERS = \ … … 68 62 # psphotApplyPSF.c 69 63 # psphotFitGalaxies.c 64 # psphotImageStats.c 65 # psphotImageBackground.c 66 # psphotOutput.c 67 # psphotLoadPixels.c -
trunk/psphot/src/psModulesUtils.c
r6522 r6571 1 1 # include "psModulesUtils.h" 2 3 // extract config informatin from config data or from image header, if specified4 // XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme5 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {6 7 char *source;8 char *keyword;9 psF32 value;10 psMetadataItem *item;11 12 // look for the entry in the config collection13 item = psMetadataLookup (config, name);14 if (item == NULL) {15 psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);16 value = psMetadataLookupF32 (status, header, name);17 if (*status == false) {18 psAbort ("pmConfigLookupF32", "no key %s in header", name);19 }20 *status = true;21 return (value);22 }23 24 // I'm either expecting a string, with the name "HD:keyword"...25 if (item->type == PS_DATA_STRING) {26 source = item->data.V;27 if (!strncasecmp (source, "HD:", 3)) {28 keyword = &source[3];29 value = psMetadataLookupF32 (status, header, keyword);30 if (*status == false) {31 psAbort ("pmConfigLookupF32", "no key %s in config", name);32 }33 *status = true;34 // psFree (item);35 return (value);36 }37 }38 39 // ... or a value (F32?)40 if (item->type == PS_DATA_F32) {41 value = item->data.F32;42 // psFree (item);43 return (value);44 }45 46 *status = false;47 psError(PS_ERR_UNKNOWN, true, "invalid item");48 return (0);49 }50 51 bool pmModelFitStatus (pmModel *model) {52 53 bool status;54 55 pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);56 status = statusFunc (model);57 58 return (status);59 }60 61 pmModel *pmModelSelect (pmSource *source) {62 switch (source->type) {63 case PM_SOURCE_STAR:64 return source->modelPSF;65 66 case PM_SOURCE_EXTENDED:67 return source->modelEXT;68 69 default:70 return NULL;71 }72 return NULL;73 }74 75 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {76 77 float obsSum = 0;78 float fitSum = 0;79 float sky = model->params->data.F32[0];80 81 *fitMag = 99.0;82 *obsMag = 99.0;83 84 pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);85 fitSum = modelFluxFunc (model->params);86 if (fitSum <= 0) return false;87 if (!isfinite(fitSum)) return false;88 *fitMag = -2.5*log10(fitSum);89 90 for (int ix = 0; ix < image->numCols; ix++) {91 for (int iy = 0; iy < image->numRows; iy++) {92 if (mask->data.U8[iy][ix]) continue;93 obsSum += image->data.F32[iy][ix] - sky;94 }95 }96 if (obsSum <= 0) return false;97 *obsMag = -2.5*log10(obsSum);98 99 return (true);100 }101 102 // force the fitted sky to meet the source flux at outer radius103 bool pmModelSkyOffset (pmModel *model, float x, float y, float radius) {104 105 float flux;106 107 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);108 109 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);110 psVector *params = model->params;111 112 coord->data.F32[0] = x + radius;113 coord->data.F32[1] = y;114 flux = modelFunc (NULL, params, coord);115 params->data.F32[0] = flux;116 117 psFree (coord);118 119 return true;120 }121 122 pmModel *pmModelCopy (pmModel *model) {123 124 pmModel *new = pmModelAlloc (model->type);125 126 new->chisq = model->chisq;127 new->nDOF = model->nDOF;128 new->nIter = model->nIter;129 new->status = model->status;130 new->radius = model->radius;131 132 for (int i = 0; i < new->params->n; i++) {133 new->params->data.F32[i] = model->params->data.F32[i];134 new->dparams->data.F32[i] = model->dparams->data.F32[i];135 }136 137 return (new);138 }139 140 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) {141 142 int Xs, Xe, Ys, Ye;143 int xi, xj, yi, yj;144 int xIs, xJs, yIs, yJs;145 int xIe, yIe;146 float flux, wt;147 148 psImage *Pi = Mi->pixels;149 psImage *Pj = Mj->pixels;150 151 psImage *Wi = Mi->weight;152 153 psImage *Ti = Mi->mask;154 psImage *Tj = Mj->mask;155 156 Xs = PS_MAX (Pi->col0, Pj->col0);157 Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);158 159 Ys = PS_MAX (Pi->row0, Pj->row0);160 Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);161 162 xIs = Xs - Pi->col0;163 xJs = Xs - Pj->col0;164 yIs = Ys - Pi->row0;165 yJs = Ys - Pj->row0;166 167 xIe = Xe - Pi->col0;168 yIe = Ye - Pi->row0;169 170 // note that this is addressing the same image pixels,171 // though only if both are source not model images172 flux = 0;173 for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {174 for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {175 if (Ti->data.U8[yi][xi]) continue;176 if (Tj->data.U8[yj][xj]) continue;177 wt = Wi->data.F32[yi][xi];178 if (wt > 0) {179 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;180 }181 }182 }183 return (flux);184 }185 186 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {187 188 int Xs, Xe, Ys, Ye;189 int xi, xj, yi, yj;190 int xIs, xJs, yIs, yJs;191 int xIe, yIe;192 float flux, wt;193 194 psImage *Pi = Mi->pixels;195 psImage *Pj = Mj->pixels;196 197 psImage *Wi = Mi->weight;198 199 psImage *Ti = Mi->mask;200 psImage *Tj = Mj->mask;201 202 Xs = PS_MAX (Pi->col0, Pj->col0);203 Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);204 205 Ys = PS_MAX (Pi->row0, Pj->row0);206 Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);207 208 xIs = Xs - Pi->col0;209 xJs = Xs - Pj->col0;210 yIs = Ys - Pi->row0;211 yJs = Ys - Pj->row0;212 213 xIe = Xe - Pi->col0;214 yIe = Ye - Pi->row0;215 216 // note that this is addressing the same image pixels,217 // though only if both are source not model images218 flux = 0;219 for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {220 for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {221 if (Ti->data.U8[yi][xi]) continue;222 if (Tj->data.U8[yj][xj]) continue;223 wt = Wi->data.F32[yi][xi];224 if (wt > 0) {225 flux += 1.0 / wt;226 }227 }228 }229 return (flux);230 }231 2 232 3 static void ppConfigFree (ppConfig *config) { -
trunk/psphot/src/psModulesUtils.h
r6522 r6571 23 23 psMetadata *recipe; // The recipe (i.e., specific setups) 24 24 psMetadata *arguments; // Command-line arguments 25 psMetadata *options; // Command-line recipe options26 25 psDB *database; // Database handle 27 26 } ppConfig; -
trunk/psphot/src/psphot.c
r6481 r6571 10 10 11 11 // load command-line arguments, options, and system config data 12 p pConfig *config = psphotArguments (&argc, argv);12 pmConfig *config = psphotArguments (&argc, argv); 13 13 14 14 // load input data (config and images (signal, noise, mask) 15 ppFile *input = psphotParseCamera (config); 16 17 // check on output mode,format, setup basic filename, etc 18 psphotOutputPrep (input, config); 15 psphotParseCamera (config); 19 16 20 17 // call psphot for each readout 21 psphotImageLoop ( input,config);18 psphotImageLoop (config); 22 19 23 20 psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); 24 21 25 psFree (input);26 22 psFree (config); 27 23 psphotCleanup (); -
trunk/psphot/src/psphot.h
r6481 r6571 10 10 # include <pmPSFtry.h> 11 11 # include <pmModelGroup.h> 12 # include "psLibUtils.h"13 # include "psModulesUtils.h"14 # include "psSparse.h"12 // # include "psLibUtils.h" 13 // # include "psModulesUtils.h" 14 // # include "psSparse.h" 15 15 # include "pmFPAConstruct.h" 16 16 # include "pmConcepts.h" … … 21 21 22 22 // top-level psphot functions 23 p pConfig *psphotArguments (int *argc, char **argv);24 ppFile *psphotParseCamera (ppConfig *config);25 bool psphotImageLoop (p pFile *input, ppConfig *config);23 pmConfig *psphotArguments (int *argc, char **argv); 24 bool psphotParseCamera (pmConfig *config); 25 bool psphotImageLoop (pmConfig *config); 26 26 27 27 bool psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe); 28 28 bool psphotReadout (pmReadout *readout, psMetadata *config); 29 bool ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum);30 29 void psphotCleanup (void); 30 31 // bool ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum); 31 32 32 33 // psphotReadout functions … … 47 48 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf); 48 49 49 // XXX deprecate these?50 // bool psphotFullFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);51 // bool psphotApplyPSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);52 // bool psphotFitExtended (pmReadout *readout, psMetadata *config, psArray *sources, psStats *skyStats);53 54 50 // basic support functions 55 pmModel *pmModelCopy (pmModel *model);56 pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);57 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);58 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);59 psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);60 51 void psphotModelGroupInit (void); 61 52 int psphotSortBySN (const void **a, const void **b); … … 63 54 bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf); 64 55 void psphotTestArguments (int *argc, char **argv); 65 bool p mCellSetMask(pmCell *cell, psMetadata *recipe);56 bool psphotMaskCell (pmCell *cell, psMetadata *recipe); 66 57 bool psphotBackgroundNames (psMetadata *arguments); 67 58 bool psphotSkyReplace (pmReadout *readout, psImage *background); … … 73 64 bool psphotInitRadiusEXT (psMetadata *config, pmModelType type); 74 65 bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model); 75 bool psphotDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius);76 bool psphotRedefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius);77 66 78 67 // output functions 68 # if (0) 79 69 bool pmSourcesWriteSX (psArray *sources, char *filename); 80 70 bool pmSourcesWriteOBJ (psArray *sources, char *filename); … … 93 83 int psphotSaveImage (psMetadata *header, psImage *image, char *filename); 94 84 bool psphotUpdateHeader (psMetadata *header, psMetadata *config); 95 int pmSourcesDophotType (pmSource *source); 96 bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key); 85 # endif 97 86 98 char *psphotSplitName (psMetadata *header); 99 void psphotOutputPrep (ppFile *file, ppConfig *config); 100 void psphotOutputCleanup (void); 101 char *psphotNameSubstitute (char *input, char *replace, char *key); 87 // void psphotOutputPrep (ppFile *file, ppConfig *config); 88 // void psphotOutputCleanup (void); 102 89 103 90 bool psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf); … … 107 94 bool psphotEvalDBL (pmSource *source, pmModel *model); 108 95 bool psphotEvalEXT (pmSource *source, pmModel *model); 109 110 // functions to support simultaneous multi-source fitting111 bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);112 bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);113 psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);114 bool pmModelFitSetInit (pmModelType type);115 void pmModelFitSetClear (void);116 96 117 97 // functions to support the source fitting process … … 124 104 psArray *psphotFitDBL (pmReadout *readout, pmSource *source); 125 105 126 // bicubic interpolation 127 psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); 128 psPlane psImageBicubeMin (psPolynomial2D *poly); 129 130 bool psImageJpegColormap (char *name); 131 bool psImageJpeg (psImage *image, char *filename, float zero, float scale); 106 // XXX these can probably be dropped: 132 107 133 108 // optional mode for clip fit? … … 143 118 const psVector *z, 144 119 const psVector *t); 120 121 // functions to support simultaneous multi-source fitting 122 bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); 123 bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); 124 psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); 125 bool pmModelFitSetInit (pmModelType type); 126 void pmModelFitSetClear (void); 127 128 // bicubic interpolation 129 psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); 130 psPlane psImageBicubeMin (psPolynomial2D *poly); 131 132 bool psImageJpegColormap (char *name); 133 bool psImageJpeg (psImage *image, char *filename, float zero, float scale); 134 135 136 // XXX deprecate 137 138 // char *psphotNameSubstitute (char *input, char *replace, char *key); 139 // char *psphotSplitName (psMetadata *header); 140 // pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); 141 // float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); 142 // float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj); 143 // psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); -
trunk/psphot/src/psphotArguments.c
r6522 r6571 1 1 # include "psphot.h" 2 # include <glob.h> 2 3 3 4 static void usage (void); 4 5 5 p pConfig *psphotArguments (int *argc, char **argv) {6 pmConfig *psphotArguments (int *argc, char **argv) { 6 7 7 8 int N; … … 12 13 psLogSetFormat ("M"); 13 14 14 ppConfig *config = ppConfigAlloc (); 15 // these other options override the PSPHOT recipe options 16 psMetadata *options = psMetadataAlloc (); 15 17 16 // load config data from default locations 17 if (!pmConfigRead(&config->site, &config->camera, &config->recipe, argc, argv, PSPHOT_RECIPE)) { 18 psErrorStackPrint(stderr, "Can't find site configuration!\n"); 19 exit (EXIT_FAILURE); 18 // run the test model (requires X,Y coordinate) 19 if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { 20 psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT", 0, "", true); 21 psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); 22 psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); 23 24 psArgumentRemove (N, argc, argv); 25 psArgumentRemove (N, argc, argv); 26 psArgumentRemove (N, argc, argv); 27 28 // specify the modeltest model 29 if ((N = psArgumentGet (*argc, argv, "-model"))) { 30 psArgumentRemove (N, argc, argv); 31 psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]); 32 psArgumentRemove (N, argc, argv); 33 } 34 35 // specify the test fit mode 36 if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { 37 psArgumentRemove (N, argc, argv); 38 psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]); 39 psArgumentRemove (N, argc, argv); 40 } 41 if ((N = psArgumentGet (*argc, argv, "-fitset"))) { 42 psArgumentRemove (N, argc, argv); 43 psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]); 44 psArgumentRemove (N, argc, argv); 45 } 20 46 } 21 47 22 // Parse other command-line arguments 23 config->arguments = psMetadataAlloc (); 24 25 // arguments (must be supplied for each run, or not used) 26 // mask image (input) is used by psphotImageLoop 27 if ((N = psArgumentGet (*argc, argv, "-mask"))) { 48 // photcode : used in output to supplement header data (argument or recipe?) 49 if ((N = psArgumentGet (*argc, argv, "-photcode"))) { 28 50 psArgumentRemove (N, argc, argv); 29 psMetadataAddStr ( config->arguments, PS_LIST_TAIL, "MASK_IMAGE", PS_META_REPLACE, "", argv[N]);51 psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); 30 52 psArgumentRemove (N, argc, argv); 31 53 } 32 54 33 // weight image (input) is used by psphotImageLoop34 if ((N = psArgumentGet (*argc, argv, "- weight"))) {55 // break : used from recipe throughout psphotReadout 56 if ((N = psArgumentGet (*argc, argv, "-break"))) { 35 57 psArgumentRemove (N, argc, argv); 36 psMetadataAddStr ( config->arguments, PS_LIST_TAIL, "WEIGHT_IMAGE", PS_META_REPLACE, "", argv[N]);58 psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); 37 59 psArgumentRemove (N, argc, argv); 38 60 } 39 61 40 // residual image (output) is used by psphotOutput41 if ((N = psArgumentGet (*argc, argv, "- resid"))) {62 // fitmode : used from recipe throughout psphotReadout 63 if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { 42 64 psArgumentRemove (N, argc, argv); 43 psMetadataAddStr ( config->arguments, PS_LIST_TAIL, "RESID_IMAGE", PS_META_REPLACE, "", argv[N]);65 psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]); 44 66 psArgumentRemove (N, argc, argv); 45 67 } 46 68 47 // background image (output) is used by psphotOutput48 if ((N = psArgumentGet (*argc, argv, "- background"))) {69 // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF 70 if ((N = psArgumentGet (*argc, argv, "-region"))) { 49 71 psArgumentRemove (N, argc, argv); 50 psMetadataAddStr ( config->arguments, PS_LIST_TAIL, "BACKGROUND_IMAGE", PS_META_REPLACE, "", argv[N]);72 psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); 51 73 psArgumentRemove (N, argc, argv); 52 74 } 53 75 54 // background model (output) is used by psphotOutput55 if ((N = psArgumentGet (*argc, argv, "-backmodel"))) {76 // other arbitrary recipe values: -D key value (all added as string) 77 while ((N = psArgumentGet (*argc, argv, "-D"))) { 56 78 psArgumentRemove (N, argc, argv); 57 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_MODEL", PS_META_REPLACE, "", argv[N]); 79 psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); 80 psArgumentRemove (N, argc, argv); 58 81 psArgumentRemove (N, argc, argv); 59 82 } 83 84 // other arbitrary recipe values: -Df key value (all added as float) 85 while ((N = psArgumentGet (*argc, argv, "-Df"))) { 86 psArgumentRemove (N, argc, argv); 87 psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); 88 psArgumentRemove (N, argc, argv); 89 psArgumentRemove (N, argc, argv); 90 } 91 92 // other arbitrary recipe values: -Di key value (all added as int) 93 while ((N = psArgumentGet (*argc, argv, "-Di"))) { 94 psArgumentRemove (N, argc, argv); 95 psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); 96 psArgumentRemove (N, argc, argv); 97 psArgumentRemove (N, argc, argv); 98 } 99 100 // load config data from default locations 101 pmConfig *config = pmConfigRead(argc, argv, PSPHOT_RECIPE); 102 103 // Storage for other command-line arguments 104 config->arguments = psMetadataAlloc (); 105 106 // save these recipe options until we have loaded the options 107 psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "PSPHOT.OPTIONS", PS_META_REPLACE, "", options); 60 108 61 109 // chip selection is used to limit chips to be processed … … 66 114 } 67 115 68 // run the test model (requires X,Y coordinate) 69 if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { 70 psMetadataAddBool (config->arguments, PS_LIST_TAIL, "TEST_FIT", 0, "", true); 71 psMetadataAddF32 (config->arguments, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); 72 psMetadataAddF32 (config->arguments, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); 116 // we load all input files onto a psArray, to be parsed later 117 psArray *input = psArrayAlloc (16); 73 118 119 // load the list of filenames the supplied file (may be a glob: "file*.fits") 120 if ((N = psArgumentGet (*argc, argv, "-file"))) { 121 glob_t globList; 74 122 psArgumentRemove (N, argc, argv); 75 psArgumentRemove (N, argc, argv); 76 psArgumentRemove (N, argc, argv); 77 78 // specify the modeltest model 79 if ((N = psArgumentGet (*argc, argv, "-model"))) { 80 psArgumentRemove (N, argc, argv); 81 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]); 82 psArgumentRemove (N, argc, argv); 123 globList.gl_offs = 0; 124 glob (argv[N], 0, NULL, &globList); 125 for (int i = 0; i < globList.gl_pathc; i++) { 126 char *filename = psStringCopy (globList.gl_pathv[i]); 127 psArrayAdd (input, 16, filename); 83 128 } 84 85 // specify the test fit mode86 if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {87 psArgumentRemove (N, argc, argv);88 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]);89 psArgumentRemove (N, argc, argv);90 }91 if ((N = psArgumentGet (*argc, argv, "-fitset"))) {92 psArgumentRemove (N, argc, argv);93 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]);94 psArgumentRemove (N, argc, argv);95 }96 }97 98 // these other options override the recipe options99 config->options = psMetadataAlloc ();100 101 // PSF : optional psf file to be loaded102 if ((N = psArgumentGet (*argc, argv, "-psf"))) {103 psArgumentRemove (N, argc, argv);104 psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_INPUT_FILE", PS_META_REPLACE, "", argv[N]);105 129 psArgumentRemove (N, argc, argv); 106 130 } 107 131 108 // PSF : optional psf file to be loaded 109 if ((N = psArgumentGet (*argc, argv, "-psfout"))) { 132 // load the list from the supplied text file 133 if ((N = psArgumentGet (*argc, argv, "-list"))) { 134 int nItems; 135 char line[1024]; // XXX limits the list lines to 1024 chars 136 char word[1024]; 137 char *filename; 138 110 139 psArgumentRemove (N, argc, argv); 111 psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_OUTPUT_FILE", PS_META_REPLACE, "", argv[N]); 140 FILE *f = fopen (argv[N], "r"); 141 if (f == NULL) { 142 psAbort ("psphot", "unable to open specified list file"); 143 } 144 while (fgets (line, 1024, f) != NULL) { 145 nItems = sscanf (line, "%s", word); 146 switch (nItems) { 147 case 0: 148 break; 149 case 1: 150 filename = psStringCopy (word); 151 psArrayAdd (input, 16, filename); 152 break; 153 default: 154 // rigid format, no comments allowed? 155 psAbort ("psphot", "error parsing input list file"); 156 break; 157 } 158 } 112 159 psArgumentRemove (N, argc, argv); 113 } 160 } 161 if (input->n < 1) usage (); 114 162 115 // photcode : used in output to supplement header data (argument or recipe?) 116 if ((N = psArgumentGet (*argc, argv, "-photcode"))) { 117 psArgumentRemove (N, argc, argv); 118 psMetadataAddStr (config->options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); 119 psArgumentRemove (N, argc, argv); 120 } 163 // input list gets places as an array on the config->arguements list 164 psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "INPUT", PS_META_REPLACE, "", input); 121 165 122 // break : used from recipe throughout psphotReadout 123 if ((N = psArgumentGet (*argc, argv, "-break"))) { 124 psArgumentRemove (N, argc, argv); 125 psMetadataAddStr (config->options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); 126 psArgumentRemove (N, argc, argv); 127 } 128 129 // fitmode : used from recipe throughout psphotReadout 130 if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { 131 psArgumentRemove (N, argc, argv); 132 psMetadataAddStr (config->options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]); 133 psArgumentRemove (N, argc, argv); 134 } 135 136 // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF 137 if ((N = psArgumentGet (*argc, argv, "-region"))) { 138 psArgumentRemove (N, argc, argv); 139 psMetadataAddStr (config->options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); 140 psArgumentRemove (N, argc, argv); 141 } 142 143 // other arbitrary recipe values: -D key value (all added as string) 144 while ((N = psArgumentGet (*argc, argv, "-D"))) { 145 psArgumentRemove (N, argc, argv); 146 psMetadataAddStr (config->options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); 147 psArgumentRemove (N, argc, argv); 148 psArgumentRemove (N, argc, argv); 149 } 150 151 // other arbitrary recipe values: -Df key value (all added as float) 152 while ((N = psArgumentGet (*argc, argv, "-Df"))) { 153 psArgumentRemove (N, argc, argv); 154 psMetadataAddF32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); 155 psArgumentRemove (N, argc, argv); 156 psArgumentRemove (N, argc, argv); 157 } 158 159 // other arbitrary recipe values: -Di key value (all added as int) 160 while ((N = psArgumentGet (*argc, argv, "-Di"))) { 161 psArgumentRemove (N, argc, argv); 162 psMetadataAddS32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); 163 psArgumentRemove (N, argc, argv); 164 psArgumentRemove (N, argc, argv); 165 } 166 167 if (*argc != 3) usage (); 166 if (*argc != 2) usage (); 168 167 169 168 // input and output positions are fixed 170 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "INPUT_FILE", PS_META_REPLACE, "", argv[1]); 171 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT_ROOT", PS_META_REPLACE, "", argv[2]); 169 psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", PS_META_REPLACE, "", argv[1]); 172 170 173 171 psTrace(__func__, 1, "Done with psphotArguments...\n"); … … 176 174 177 175 static void usage (void) { 178 fprintf (stderr, "USAGE: psphot (image)(output)\n");176 fprintf (stderr, "USAGE: psphot [-file image(s)] [-list imagelist] (output)\n"); 179 177 exit (2); 180 178 } -
trunk/psphot/src/psphotChoosePSF.c
r6481 r6571 3 3 // 2006.02.07 : no leaks! 4 4 // try PSF models and select best option 5 pmPSF *psphotChoosePSF (ps Metadata *config, psArray *sources) {5 pmPSF *psphotChoosePSF (psArray *sources, psMetadata *config) { 6 6 7 7 bool status; -
trunk/psphot/src/psphotCleanup.c
r6399 r6571 9 9 psTimeFinalize (); 10 10 pmConceptsDone (); 11 fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");11 fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); 12 12 return; 13 13 } -
trunk/psphot/src/psphotEnsemblePSF.c
r6481 r6571 4 4 // 2006.02.07 : no leaks! 5 5 // fit all reasonable sources with the linear PSF model 6 bool psphotEnsemblePSF (pmReadout *readout, ps Metadata *config, psArray *sources, pmPSF *psf, bool final) {6 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf, bool final) { 7 7 8 8 bool status; -
trunk/psphot/src/psphotImageMedian.c
r6481 r6571 9 9 static int MAX_SAMPLE_PIXELS; 10 10 11 static char *backImage = NULL;12 static char *backModel = NULL;13 14 bool psphotBackgroundNames (psMetadata *arguments) {15 16 bool status;17 18 backImage = psMetadataLookupStr (&status, arguments, "BACKGROUND_IMAGE");19 backModel = psMetadataLookupStr (&status, arguments, "BACKGROUND_MODEL");20 21 return true;22 }23 24 11 // generate the median in NxN boxes, clipping heavily 25 12 // linear interpolation to generate full-scale model 26 psImage *psphotImageMedian (p mReadout *readout, psMetadata *config)13 psImage *psphotImageMedian (psMetadata *config, pmFPAview *view) 27 14 { 28 15 bool status; 29 16 psRegion region; 17 psImage *model = NULL; 18 19 psTimerStart ("psphot"); 20 21 // find the currently selected readout 22 psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); 23 pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); 24 pmReadout *readout = pmFPAviewThisReadout (view, input); 25 26 psImage *image = readout->image; 27 psImage *mask = readout->mask; 30 28 31 29 rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); 32 30 MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, config, "IMSTATS_NPIX"); 33 31 if (!status) MAX_SAMPLE_PIXELS = 1000; 34 35 psTimerStart ("psphot");36 37 psImage *image = readout->image;38 psImage *mask = readout->mask;39 32 40 33 // dimensions of input & output image … … 56 49 int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY; 57 50 58 // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);59 psImage *model = psImageAlloc (nx, ny, PS_TYPE_F32);51 // select model pixels, from output background model file, or create 52 model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32); 60 53 61 54 // measure clipped median for subimages … … 76 69 77 70 // XXX psImageStats is still very inefficient and poorly coded... 71 // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); 78 72 // stats = psImageStats (stats, subset, maskset, 0xff); 79 // model->data.F32[iy][ix] = stats->clippedMean;73 // backMdl->data.F32[iy][ix] = stats->clippedMean; 80 74 } 81 75 } … … 83 77 84 78 // linear interpolation to full-scale 85 86 psImage *background = psImageAlloc (Nx, Ny, PS_TYPE_F32); 79 80 // select background pixels, from output background file, or create 81 background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32); 87 82 88 83 // XXX this code skips the initial pixels … … 216 211 psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot")); 217 212 213 // back-sub image pixels, from output background file (don't create if not requested) 214 backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32); 215 218 216 // subtract the background model 219 217 for (int j = 0; j < image->numRows; j++) { … … 221 219 if (!mask->data.U8[j][i]) { 222 220 image->data.F32[j][i] -= background->data.F32[j][i]; 223 } 224 } 225 } 226 227 // optionally save background 228 if (backImage != NULL) psphotSaveImage (NULL, background, backImage); 229 if (backModel != NULL) psphotSaveImage (NULL, model, backModel); 221 if (backSub) { 222 backSub->data.F32[j][i] = image->data.F32[j][i]; 223 } 224 } 225 } 226 } 230 227 231 228 psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot")); 232 229 psFree (rnd); 230 231 // it is safe to free all of these: either they are saved on config->files or are local 233 232 psFree (model); 234 return (background); 233 psFree (background); 234 psFree (backSub); 235 return true; 235 236 } 236 237 -
trunk/psphot/src/psphotMagnitudes.c
r6495 r6571 1 1 # include "psphot.h" 2 3 // XXX EAM : the apMag should only be calculated for the brighter sources?4 // XXX EAM : SN limit set by user?5 // XXX EAM : masked region should be (optionally) elliptical6 pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) {7 8 int status;9 bool isPSF;10 float x, y;11 float rflux;12 float radius;13 pmModel *model;14 15 switch (source->type) {16 case PM_SOURCE_STAR:17 model = source->modelPSF;18 if (model == NULL) return NULL;19 radius = (apRadius > 0) ? apRadius : model->radius;20 isPSF = true;21 break;22 23 case PM_SOURCE_EXTENDED:24 model = source->modelEXT;25 if (model == NULL) return NULL;26 radius = model->radius;27 isPSF = false;28 break;29 30 default:31 return NULL;32 }33 34 x = model->params->data.F32[2];35 y = model->params->data.F32[3];36 37 // replace source flux38 pmSourceAddModel (source->pixels, source->mask, model, false, false);39 40 // set aperture mask circle of PSF_FIT_RADIUS41 psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED);42 43 // measure object photometry44 status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask);45 46 // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS47 if (isPSF && (psf != NULL)) {48 if (psf->growth != NULL) {49 source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius);50 }51 52 rflux = pow (10.0, 0.4*source->fitMag);53 source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux;54 source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);55 }56 57 // unmask aperture58 psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED);59 60 // subtract object, leave local sky61 pmSourceSubModel (source->pixels, source->mask, model, false, false);62 63 if (!status) return NULL;64 return model;65 }66 67 /*68 aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)69 (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)70 (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)71 72 */73 2 74 3 bool psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf) { … … 88 17 } 89 18 19 // XXX should this function use RADIUS or should it -
trunk/psphot/src/psphotModelGroupInit.c
r5718 r6571 1 1 # include "psphot.h" 2 3 // Add locally-defined models here. As these mature, they can be moved to 4 // psModule/src/objects/models 2 5 3 6 # include "models/pmModel_TAUSS.c" … … 6 9 {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, 7 10 }; 8 9 11 10 12 void psphotModelGroupInit (void) -
trunk/psphot/src/psphotModelTest.c
r6481 r6571 3 3 static char DEFAULT_MODE[] = "EXT"; 4 4 5 bool psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe) {5 bool psphotModelTest (pmReadout *readout, recipe) { 6 6 7 7 bool status; … … 15 15 16 16 // run model fitting tests on a single source 17 if (!psMetadataLookupBool (&status, arguments, "TEST_FIT")) return false;17 if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false; 18 18 19 19 // what fitting mode to use? 20 char *psfModeWord = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODE");20 char *psfModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE"); 21 21 if (!status) { 22 22 psfModeWord = DEFAULT_MODE; … … 33 33 } else { 34 34 // find the model: supplied by user or first in the PSF_MODEL list 35 char *modelName = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODEL");35 char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); 36 36 if (modelName == NULL) { 37 37 // get the list pointers for the PSF_MODEL entries … … 57 57 58 58 // find the fitting parameters (try test values first) 59 float INNER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_INNER_RADIUS");59 float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS"); 60 60 if (!status) { 61 61 INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); 62 62 } 63 63 64 float OUTER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_OUTER_RADIUS");64 float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS"); 65 65 if (!status) { 66 66 OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); 67 67 } 68 68 69 float RADIUS = psMetadataLookupF32 (&status, arguments, "TEST_FIT_RADIUS");69 float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS"); 70 70 if (!status) { 71 71 RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); 72 72 } 73 73 74 float mRADIUS = psMetadataLookupF32 (&status, arguments, "TEST_MOMENTS_RADIUS");74 float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS"); 75 75 if (!status) { 76 76 mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); … … 78 78 79 79 // define the source of interest 80 float xObj = psMetadataLookupF32 (&status, arguments, "TEST_FIT_X");81 float yObj = psMetadataLookupF32 (&status, arguments, "TEST_FIT_Y");80 float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); 81 float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); 82 82 83 83 // construct the source structures … … 120 120 } 121 121 sprintf (name, "TEST_FIT_PAR%d", i); 122 value = psMetadataLookupF32 (&status, arguments, name);122 value = psMetadataLookupF32 (&status, recipe, name); 123 123 if (status) { 124 124 params[i] = value; … … 154 154 psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED); 155 155 156 char *fitset = psMetadataLookupStr (&status, arguments, "TEST_FIT_SET");156 char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); 157 157 if (status) { 158 158 status = psphotFitSet (source, model, fitset, psfMode); -
trunk/psphot/src/psphotOutput.c
r6522 r6571 1 1 # include "psphot.h" 2 3 static int extNumber = -1;4 5 static char *outputName = NULL;6 static char *outputRoot = NULL;7 static char *outputMode = NULL;8 static char *outputFormat = NULL;9 10 static char *extNumberFormat = NULL;11 static char *extNameKey = NULL;12 static char *extRoot = NULL;13 static char *psfFile = NULL;14 static char *psfSample = NULL;15 16 static psFits *outputFits = NULL;17 static FILE *outputFile = NULL;18 19 bool psphotDataIO (psphotData *data, ppConfig *config, int chip, int cell) {20 21 // load input data, if depth is appropriate22 ppImageLoadDepth loadDepth = ppImageCheckDepth (config->recipe, "LOAD.DEPTH");23 if ((chip == -1) && (cell == -1) && (loadDepth == PP_LOAD_FPA)) {24 psTrace(__func__, 1, "Loading pixels for FPA...\n");25 ppImageLoadPixels(data->input, config->database, -1, -1);26 }27 if ((chip != -1) && (cell == -1) && (loadDepth == PP_LOAD_CHIP)) {28 psTrace(__func__, 1, "Loading pixels for chip %d...\n", i);29 ppImageLoadPixels(data->input, config->database, chip, -1);30 }31 if ((chip != -1) && (cell != -1) && (loadDepth == PP_LOAD_CELL)) {32 psTrace(__func__, 1, "Loading pixels for chip %d, cell %d...\n", i, j);33 ppImageLoadPixels(data->input, config->database, chip, cell);34 }35 36 37 38 if (saveDepth == PP_LOAD_FPA) {39 psTrace(__func__, 1, "Saving objects for FPA...\n");40 psphotOutputOpen (file->fpa->phu);41 }42 43 ppImageLoadDepth saveDepth = ppImageCheckDepth (config->recipe, "SAVE.DEPTH");44 45 46 }47 48 void psphotOutputPrep (ppFile *file, ppConfig *config) {49 50 bool status;51 52 outputRoot = psMetadataLookupStr (&status, config->arguments, "OUTPUT_ROOT");53 if (!status) psAbort ("psphot", "output file not specified");54 55 outputName = psMetadataLookupStr (&status, config->recipe, "OUTPUT_NAME");56 outputMode = psMetadataLookupStr (&status, config->recipe, "OUTPUT_MODE");57 outputFormat = psMetadataLookupStr (&status, config->recipe, "OUTPUT_FORMAT");58 59 psfFile = psMetadataLookupStr (&status, config->recipe, "PSF_OUTPUT_FILE");60 psfSample = psMetadataLookupStr (&status, config->recipe, "PSF_SAMPLE_FILE");61 62 // rules to construct output names63 extNumberFormat = psMetadataLookupStr (&status, config->recipe, "EXTNUMBER_FORMAT");64 extNameKey = psMetadataLookupStr (&status, config->recipe, "EXTNAME");65 extRoot = psMetadataLookupStr (&status, config->recipe, "EXTROOT");66 67 headExtname = psMetadataLookupStr (&status, config->recipe, "HEAD_EXTNAME");68 dataExtname = psMetadataLookupStr (&status, config->recipe, "DATA_EXTNAME");69 70 if (extNumberFormat == NULL) {71 extNumberFormat = psStringCopy ("%02d");72 }73 74 psStringStrip (outputName);75 psStringStrip (extNameKey);76 psStringStrip (extRoot);77 78 // validate the outputMode and outputFormat79 status = false;80 status |= !strcasecmp (outputFormat, "TEXT");81 status |= !strcasecmp (outputFormat, "OBJ");82 status |= !strcasecmp (outputFormat, "SX");83 status |= !strcasecmp (outputFormat, "CMP");84 status |= !strcasecmp (outputFormat, "CMF");85 if (!status) {86 psLogMsg ("psphot", PS_LOG_ERROR, "invalid output format %s", outputFormat);87 exit (1);88 }89 90 // validate the outputMode and outputFormat91 status = false;92 status |= !strcasecmp (outputMode, "MEF");93 status |= !strcasecmp (outputMode, "SPLIT");94 if (!status) {95 psLogMsg ("psphot", PS_LOG_ERROR, "invalid output mode %s", outputMode);96 exit (1);97 }98 99 // for MEF output, we only allow certain output formats100 if (!strcasecmp (outputMode, "MEF")) {101 if (!strcasecmp (outputFormat, "CMF")) goto valid_format;102 psLogMsg ("psphot", PS_LOG_ERROR, "MEF output mode unavailable");103 exit (1);104 }105 valid_format:106 107 // register background image-file names108 psphotBackgroundNames (config->arguments);109 110 // save so freeing config will not drop these references111 psMemCopy (outputRoot);112 psMemCopy (outputName);113 psMemCopy (outputMode);114 psMemCopy (outputFormat);115 psMemCopy (extNumberFormat);116 psMemCopy (extNameKey);117 psMemCopy (extRoot);118 psMemCopy (psfFile);119 psMemCopy (psfSample);120 121 return;122 }123 124 bool psphotOutputOpen (psMetadata *phu){125 126 /* outputFile is used by text-format outputs */127 128 psFree (outputFile);129 outputFile = psphotOutputName (phu, outputName);130 131 /* for fits-format output, open outputFits */132 if (!strcasecmp (outputFormat, "CMF")) {133 if (outputFits != NULL) {134 psFitsClose (outputFits);135 }136 outputFits = psFitsOpen (outputFile, "w");137 138 if (!strcasecmp (outputMode, "MEF")) {139 psFitsWriteHeader (outputFits, phu);140 }141 }142 return true;143 }144 145 bool psphotOutputClose (void) {146 147 if (outputFile != NULL) {148 psFitsClose (outputFile);149 }150 return true;151 }152 2 153 3 // output functions: we have several fixed modes (sx, obj, cmp) … … 159 9 psMetadata *header = pmReadoutGetHeader (readout); 160 10 161 psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");162 11 pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); 163 164 // starting value is -1 165 extNumber ++; 166 167 // for SPLIT output, we need to open a file for each output call 168 if (!strcasecmp (outputMode, "SPLIT")) { 169 outputFile = psphotSplitName (header); 170 } 171 172 if (!strcasecmp (outputMode, "MEF")) { 173 headExt = psphotOutputName (header, headExtname); 174 dataExt = psphotOutputName (header, dataExtname); 175 } 176 if (!strcasecmp (outputMode, "SPLIT")) { 177 dataExt = psphotOutputName (header, dataExtname); 178 } 179 180 char *residImage = psMetadataLookupStr (&status, arguments, "RESID_IMAGE"); 12 // sample PSF images?? 181 13 if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); 182 if ( residImage != NULL) psphotSaveImage (header, readout->image, residImage);14 if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); 183 15 184 16 if (psfFile != NULL) { … … 187 19 psFree (psfData); 188 20 } 189 if (outputFile == NULL) {190 psLogMsg ("output", 3, "no data output file selected");191 return;192 }193 if (outputFormat == NULL) {194 psLogMsg ("output", 3, "no data output format selected");195 psFree (outputFile);196 return;197 }198 if (!strcasecmp (outputFormat, "SX")) {199 pmSourcesWriteSX (sources, outputFile);200 psFree (outputFile);201 return;202 }203 if (!strcasecmp (outputFormat, "OBJ")) {204 pmSourcesWriteOBJ (sources, outputFile);205 psFree (outputFile);206 return;207 }208 if (!strcasecmp (outputFormat, "CMP")) {209 pmSourcesWriteCMP (sources, outputFile, header);210 psFree (outputFile);211 return;212 }213 if (!strcasecmp (outputFormat, "CMF")) {214 pmSourcesWriteCMF (sources, header, headExt, dataExt);215 psFree (outputFile);216 return;217 }218 if (!strcasecmp (outputFormat, "TEXT")) {219 pmSourcesWriteText (sources, outputFile);220 psFree (outputFile);221 return;222 }223 224 psAbort ("psphot", "unknown output mode %s", outputFormat);225 }226 227 // dophot-style output list with fixed line width228 bool pmSourcesWriteOBJ (psArray *sources, char *filename) {229 230 int type;231 psF32 *PAR, *dPAR;232 float dmag, apResid;233 234 psTimerStart ("string");235 236 psLine *line = psLineAlloc (104); // 104 is dophot-defined line length237 238 FILE *f = fopen (filename, "w");239 if (f == NULL) {240 psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);241 return false;242 }243 244 // write sources with models245 for (int i = 0; i < sources->n; i++) {246 pmSource *source = (pmSource *) sources->data[i];247 pmModel *model = pmModelSelect (source);248 if (model == NULL) continue;249 250 PAR = model->params->data.F32;251 dPAR = model->dparams->data.F32;252 253 dmag = dPAR[1] / PAR[1];254 type = pmSourceDophotType (source);255 apResid = source->apMag - source->fitMag;256 257 psLineInit (line);258 psLineAdd (line, "%3d", type);259 psLineAdd (line, "%8.2f", PAR[2]);260 psLineAdd (line, "%8.2f", PAR[3]);261 psLineAdd (line, "%8.3f", source->fitMag);262 psLineAdd (line, "%6.3f", dmag);263 psLineAdd (line, "%9.2f", PAR[0]);264 psLineAdd (line, "%9.3f", PAR[4]);265 psLineAdd (line, "%9.3f", PAR[5]);266 psLineAdd (line, "%7.2f", PAR[6]);267 psLineAdd (line, "%8.3f", 99.999);268 psLineAdd (line, "%8.3f", source->apMag);269 psLineAdd (line, "%8.2f\n", apResid);270 fwrite (line->line, 1, line->Nline, f);271 }272 fclose (f);273 psFree (line);274 fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n);275 276 psTimerStart ("string");277 278 f = fopen ("test.obj", "w");279 if (f == NULL) {280 psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj");281 return false;282 }283 284 char *string;285 // write sources with models286 for (int i = 0; i < sources->n; i++) {287 pmSource *source = (pmSource *) sources->data[i];288 pmModel *model = pmModelSelect (source);289 if (model == NULL) continue;290 291 PAR = model->params->data.F32;292 dPAR = model->dparams->data.F32;293 294 dmag = dPAR[1] / PAR[1];295 type = pmSourceDophotType (source);296 apResid = source->apMag - source->fitMag;297 298 string = NULL;299 psStringAppend (&string, "%3d", type);300 psStringAppend (&string, "%8.2f", PAR[2]);301 psStringAppend (&string, "%8.2f", PAR[3]);302 psStringAppend (&string, "%8.3f", source->fitMag);303 psStringAppend (&string, "%6.3f", dmag);304 psStringAppend (&string, "%9.2f", PAR[0]);305 psStringAppend (&string, "%9.3f", PAR[4]);306 psStringAppend (&string, "%9.3f", PAR[5]);307 psStringAppend (&string, "%7.2f", PAR[6]);308 psStringAppend (&string, "%8.3f", 99.999);309 psStringAppend (&string, "%8.3f", source->apMag);310 psStringAppend (&string, "%8.2f\n", apResid);311 fwrite (string, 1, strlen(string), f);312 psFree (string);313 }314 fclose (f);315 fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n);316 317 return true;318 }319 320 // elixir-mode / sextractor-style output list with fixed line width321 bool pmSourcesWriteSX (psArray *sources, char *filename) {322 323 psF32 *PAR, *dPAR;324 float dmag;325 326 psLine *line = psLineAlloc (110); // 110 is sextractor line length327 328 FILE *f = fopen (filename, "w");329 if (f == NULL) {330 psLogMsg ("WriteSourceSX", 3, "can't open output file for output %s\n", filename);331 return false;332 }333 334 // write sources with models335 for (int i = 0; i < sources->n; i++) {336 pmSource *source = (pmSource *) sources->data[i];337 pmModel *model = pmModelSelect (source);338 if (model == NULL) continue;339 340 PAR = model->params->data.F32;341 dPAR = model->dparams->data.F32;342 343 // pmSourceSextractType (source, &type, &flags); XXX EAM : implement me...344 345 dmag = dPAR[1] / PAR[1];346 347 psLineInit (line);348 psLineAdd (line, "%5.2f", 0.0); // should be type349 psLineAdd (line, "%11.3f", PAR[2]);350 psLineAdd (line, "%11.3f", PAR[3]);351 psLineAdd (line, "%9.4f", source->fitMag);352 psLineAdd (line, "%9.4f", dmag);353 psLineAdd (line, "%13.4f", PAR[0]);354 psLineAdd (line, "%9.2f", 0.0); // should be FWHMx355 psLineAdd (line, "%9.2f", 0.0); // should be FWHMy356 psLineAdd (line, "%6.1f", 0.0); // should be Theta357 psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO358 psLineAdd (line, "%9.4f", source->apMag);359 psLineAdd (line, "%4d\n", 0); // should be flags360 fwrite (line->line, 1, line->Nline, f);361 }362 fclose (f);363 return true;364 }365 366 // elixir-style pseudo FITS table (header + ascii list)367 bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header) {368 369 int i, type;370 psMetadataItem *mdi;371 psF32 *PAR, *dPAR;372 float dmag, lsky;373 bool status;374 375 // find config information for output header376 float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT");377 if (!status) ZERO_POINT = 25.0;378 379 // write necessary information to output header380 psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n);381 382 // create file, write-out header383 psFits *fits = psFitsOpen (filename, "w");384 385 // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0386 mdi = psMetadataLookup (header, "NAXIS");387 mdi->data.S32 = 0;388 mdi->type = PS_DATA_S32;389 390 // psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "", 0);391 psFitsWriteHeader (header, fits);392 psFitsClose (fits);393 394 // re-open, add data to end of file395 FILE *f = fopen (filename, "a+");396 if (f == NULL) {397 psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);398 return false;399 }400 fseek (f, 0, SEEK_END);401 402 psLine *line = psLineAlloc (67); // 66 is imclean-defined line length403 404 // write sources with models first405 for (i = 0; i < sources->n; i++) {406 pmSource *source = (pmSource *) sources->data[i];407 pmModel *model = pmModelSelect (source);408 if (model == NULL) continue;409 410 PAR = model->params->data.F32;411 dPAR = model->dparams->data.F32;412 413 dmag = dPAR[1] / PAR[1];414 type = pmSourceDophotType (source);415 lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);416 417 psLineInit (line);418 psLineAdd (line, "%6.1f ", PAR[2]);419 psLineAdd (line, "%6.1f ", PAR[3]);420 psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT));421 psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag)));422 psLineAdd (line, "%2d ", type);423 psLineAdd (line, "%3.1f ", lsky);424 psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal425 psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT));426 psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x'427 psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y'428 psLineAdd (line, "%5.1f\n", 0); // should be theta429 fwrite (line->line, 1, line->Nline, f);430 }431 fclose (f);432 return true;433 }434 435 // elixir-style FITS table output (header + table in 1st extension)436 // this format consists of a header derived from the image header437 // followed by a zero-size matrix, followed by the table data438 bool pmSourcesWriteCMF (psArray *sources, char *filename, psMetadata *header) {439 440 psArray *table;441 psMetadataItem *mdi;442 psMetadata *row;443 psMetadata *theader;444 int i, type;445 psF32 *PAR, *dPAR;446 float dmag, lsky;447 bool status;448 449 // find config information for output header450 float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT");451 452 table = psArrayAlloc (sources->n);453 table->n = 0;454 455 for (i = 0; i < sources->n; i++) {456 pmSource *source = (pmSource *) sources->data[i];457 pmModel *model = pmModelSelect (source);458 if (model == NULL) continue;459 460 PAR = model->params->data.F32;461 dPAR = model->dparams->data.F32;462 463 dmag = dPAR[1] / PAR[1];464 type = pmSourceDophotType (source);465 lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);466 467 row = psMetadataAlloc ();468 psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]);469 psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]);470 psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT);471 psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag));472 psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", 32.0);473 psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", source->apMag + ZERO_POINT);474 psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky);475 psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", type);476 psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]);477 psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]);478 psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0");479 psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123");480 481 psArrayAdd (table, 100, row);482 psFree (row);483 }484 485 // write the correct number of sources in the header486 psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n);487 psMetadataAdd (header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true);488 489 // set NAXIS to 0 (we don't write out the data array)490 mdi = psMetadataLookup (header, "NAXIS");491 mdi->data.S32 = 0;492 mdi->type = PS_DATA_S32;493 494 // create the basic table header495 theader = psMetadataAlloc ();496 psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");497 498 // write out the header and table499 psFits *fits = psFitsOpen (filename, "w");500 psFitsWriteHeader (header, fits);501 psFitsWriteTable (fits, theader, table);502 psFitsClose (fits);503 504 return true;505 }506 507 /***** Text Output Methods *****/508 509 bool pmSourcesWriteText (psArray *sources, char *filename) {510 511 char *name = (char *) psAlloc (strlen(filename) + 10);512 513 sprintf (name, "%s.psf.dat", filename);514 pmModelWritePSFs (sources, name);515 516 sprintf (name, "%s.ext.dat", filename);517 pmModelWriteEXTs (sources, name);518 519 sprintf (name, "%s.nul.dat", filename);520 pmModelWriteNULLs (sources, name);521 522 sprintf (name, "%s.mnt.dat", filename);523 pmMomentsWriteText (sources, name);524 525 psFree (name);526 return true;527 }528 529 // write the PSF sources to an output file530 bool pmModelWritePSFs (psArray *sources, char *filename) {531 532 double dPos, dMag;533 int i, j;534 FILE *f;535 psF32 *PAR, *dPAR;536 pmModel *model;537 538 f = fopen (filename, "w");539 if (f == NULL) {540 psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename);541 return false;542 }543 544 // write sources with models first545 for (i = 0; i < sources->n; i++) {546 pmSource *source = (pmSource *) sources->data[i];547 if (source->type != PM_SOURCE_STAR) continue;548 model = source->modelPSF;549 if (model == NULL) continue;550 551 PAR = model->params->data.F32;552 dPAR = model->dparams->data.F32;553 554 // dPos is positional error, dMag is mag error555 dPos = hypot (dPAR[2], dPAR[3]);556 dMag = dPAR[1] / PAR[1];557 558 fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ",559 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);560 561 for (j = 4; j < model->params->n; j++) {562 fprintf (f, "%9.6f ", PAR[j]);563 }564 fprintf (f, " : ");565 for (j = 4; j < model->params->n; j++) {566 fprintf (f, "%9.6f ", dPAR[j]);567 }568 fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",569 source[0].apMag, source[0].type, source[0].mode,570 log10(model[0].chisq/model[0].nDOF),571 source[0].moments->SN,572 model[0].radius,573 model[0].nDOF,574 model[0].nIter);575 }576 fclose (f);577 return true;578 }579 580 // dump the sources to an output file581 bool pmModelWriteEXTs (psArray *sources, char *filename) {582 583 double dPos, dMag;584 int i, j;585 FILE *f;586 psF32 *PAR, *dPAR;587 pmModel *model;588 589 f = fopen (filename, "w");590 if (f == NULL) {591 psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename);592 return false;593 }594 595 // write sources with models first596 for (i = 0; i < sources->n; i++) {597 pmSource *source = (pmSource *) sources->data[i];598 599 if (source->type != PM_SOURCE_EXTENDED) continue;600 model = pmSourceMagnitudes (source, NULL, 0.0);601 if (model == NULL) continue;602 603 PAR = model->params->data.F32;604 dPAR = model->dparams->data.F32;605 606 // dPos is shape error607 // XXX these are hardwired for SGAUSS608 dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5]));609 dMag = dPAR[1] / PAR[1];610 611 fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ",612 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);613 614 for (j = 4; j < model->params->n; j++) {615 fprintf (f, "%9.6f ", PAR[j]);616 }617 fprintf (f, " : ");618 for (j = 4; j < model->params->n; j++) {619 fprintf (f, "%9.6f ", dPAR[j]);620 }621 fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",622 source->apMag,623 source[0].type, source[0].mode,624 log10(model[0].chisq/model[0].nDOF),625 source[0].moments->SN,626 model[0].radius,627 model[0].nDOF,628 model[0].nIter);629 }630 fclose (f);631 return true;632 }633 634 // dump the sources to an output file635 bool pmModelWriteNULLs (psArray *sources, char *filename)636 {637 638 int i;639 FILE *f;640 pmMoments *moment = NULL;641 pmSource *source = NULL;642 643 f = fopen (filename, "w");644 if (f == NULL) {645 psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);646 return false;647 }648 649 pmMoments *empty = pmMomentsAlloc ();650 651 // write sources with models first652 for (i = 0; i < sources->n; i++) {653 source = sources->data[i];654 655 // skip these sources (in PSF or EXT)656 if (source->type == PM_SOURCE_STAR) continue;657 if (source->type == PM_SOURCE_EXTENDED) continue;658 659 if (source->moments == NULL) {660 moment = empty;661 } else {662 moment = source->moments;663 }664 665 fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n",666 source->peak->x, source->peak->y, source->peak->counts,667 source->moments->x, source->moments->y,668 source->moments->Sx, source->moments->Sy,669 source->moments->Sum, source->moments->Peak,670 source->moments->Sky, source->moments->SN,671 source->moments->nPixels, source->type);672 }673 fclose (f);674 psFree (empty);675 return true;676 }677 678 // write the moments to an output file679 bool pmMomentsWriteText (psArray *sources, char *filename)680 {681 682 int i;683 FILE *f;684 pmSource *source = NULL;685 686 f = fopen (filename, "w");687 if (f == NULL) {688 psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename);689 return false;690 }691 692 for (i = 0; i < sources->n; i++) {693 source = sources->data[i];694 if (source->moments == NULL)695 continue;696 fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n",697 source->peak->x, source->peak->y, source->peak->counts,698 source->moments->x, source->moments->y,699 source->moments->Sx, source->moments->Sy,700 source->moments->Sum, source->moments->Peak,701 source->moments->Sky, source->moments->SN,702 source->moments->nPixels, source->type, source->mode);703 }704 fclose (f);705 return true;706 }707 708 // write the peaks to an output file709 bool pmPeaksWriteText (psArray *peaks, char *filename) {710 711 int i;712 FILE *f;713 714 f = fopen (filename, "w");715 if (f == NULL) {716 psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename);717 return false;718 }719 720 for (i = 0; i < peaks->n; i++) {721 pmPeak *peak = peaks->data[i];722 if (peak == NULL) continue;723 fprintf (f, "%5d %5d %7.1f\n",724 peak->x, peak->y, peak->counts);725 }726 fclose (f);727 return true;728 }729 730 // translations between psphot object types and dophot object types731 int pmSourceDophotType (pmSource *source) {732 733 switch (source->type) {734 735 case PM_SOURCE_DEFECT:736 case PM_SOURCE_SATURATED:737 return (8);738 739 case PM_SOURCE_STAR:740 if (source->mode & PM_SOURCE_SATSTAR) return (10);741 if (source->mode & PM_SOURCE_POOR) return (7);742 if (source->mode & PM_SOURCE_FAIL) return (4);743 return (1);744 745 case PM_SOURCE_EXTENDED:746 return (2);747 748 default:749 return (0);750 }751 return (0);752 21 } 753 22 -
trunk/psphot/src/psphotParseCamera.c
r6379 r6571 8 8 } } 9 9 10 // 2006.02.07 : no leaks! 11 ppFile *psphotParseCamera (ppConfig *config) { 10 bool psphotParseCamera (pmConfig *config) { 12 11 13 ppFile *input = ppFileAlloc (); 12 // psphot is supplied with a list of input images (may be only one image) 13 psArray *infiles = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); 14 if (infiles->n < 1) psAbort (__func__, "empty input list"); 14 15 15 input->filename = psMetadataLookupStr(NULL, config->arguments, "INPUT_FILE"); 16 psMemCopy (input->filename); // keep for external use 17 18 // Open the input image 19 psLogMsg("psphot", PS_LOG_INFO, "Opening input image: %s\n", input->filename); 20 input->fits = psFitsOpen (input->filename, "r"); // File handle for FITS file 21 if (! input->fits) { 22 psErrorStackPrint(stderr, "Can't open input image: %s\n", input->filename); 23 exit(EXIT_FAILURE); 16 // if no camera has been specified, use the first image as a template for the rest. 17 if (config->camera == NULL) { 18 psFits *fits = psFitsOpen (infiles->data[0], "r"); 19 psMetadata *phu = psFitsReadHeader (NULL, fits); 20 config->camera = pmConfigCameraFromHeader (config->site, phu); 21 psFitsClose (fits); 22 psFree (phu); 24 23 } 25 input->phu = psFitsReadHeader(NULL, input->fits); // FITS header 26 27 // Get camera configuration from header if not already defined 24 // There's no point in continuing if we can't recognize what we've got 28 25 if (! config->camera) { 29 config->camera = pmConfigCameraFromHeader(config->site, input->phu); 30 if (! config->camera) { 31 // There's no point in continuing if we can't recognize what we've got 32 psErrorStackPrint(stderr, "Can't find camera configuration!\n"); 33 exit(EXIT_FAILURE); 34 } 35 } else if (! pmConfigValidateCamera(config->camera, input->phu)) { 36 // There's no point in continuing if what we've got doesn't match what we've been told 37 psError(PS_ERR_IO, true, "%s does not seem to be from the specified camera.\n", 38 input->filename); 39 exit(EXIT_FAILURE); 26 psErrorStackPrint(stderr, "Can't find camera configuration!\n"); 27 exit(EXIT_FAILURE); 40 28 } 41 29 42 // Determine the correct recipe to use (from camera or from user) 43 // if config->recipe is not NULL, it is a user-supplied recipe 44 if (!config->recipe) { 45 config->recipe = pmConfigRecipeFromCamera(config->camera, PSPHOT_RECIPE); 46 if (config->recipe == NULL) { 47 psErrorStackPrint(stderr, "Can't find recipe configuration!\n"); 48 exit(EXIT_FAILURE); 49 } 30 // XXX place the "data" element on the pmConfig structure? 31 config->files = psMetadataAlloc (); 32 33 // build the template fpa, set up the basic view 34 pmFPA *input = pmFPAConstruct (camera); 35 36 // assign the I/O files (potentially) needed by psphot 37 pmFPAfile *file = pmFPAfileDefine (config->files, camera, input, "PSPHOT.INPUT"); 38 39 pmFPAview *view = pmFPAviewAlloc (input, config->camera, 0); 40 for (int i = 0; i < infiles->n; i++) { 41 // XXX save the phu from above? 42 psFits *fits = psFitsOpen (infiles->data[i], "r"); 43 psMetadata *phu = psFitsReadHeader (NULL, fits); 44 pmConfigValidateCamera (config->camera, phu); 45 46 // set the view to the corresponding entry for this phu 47 // XXX should I save the phu data at the appropriate location in the fpa? 48 pmFPAsetView (view, phu); 49 50 // XXX is this the correct psMD to save the filename? 51 name = pmFPAviewNameFromRule (file->filextra, view); 52 psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]); 50 53 } 54 pmFPAfileDefine (config->files, camera, input, "PSPHOT.OUTPUT"); 55 pmFPAfileDefine (config->files, camera, input, "PSPHOT.RESID"); 56 pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_INPUT"); 57 pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_OUTPUT"); 58 59 // build the template fpa, set up the basic view 60 // supply the backgnd with a different camera? 61 pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKSUB"); 62 pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKGND"); 63 pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKMDL"); 64 pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.PSF_SAMPLE"); 51 65 52 66 // recipe override values (command-line options): 53 psMetadataItem *item = NULL; 54 psMetadataIterator *iter = psMetadataIteratorAlloc (config->options, PS_LIST_HEAD, NULL); 55 while ((item = psMetadataGetAndIncrement (iter)) != NULL) { 56 psMetadataAddItem (config->recipe, item, PS_LIST_TAIL, PS_META_REPLACE); 67 psMetadata *options = psMetadataLookupPtr (status, config->arguments, "PSPHOT.OPTIONS"); 68 psMetadata *recipe = psMetadataLookupPtr (status, config->recipes, "PSPHOT"); 69 psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL); 70 while ((psMetadataItem *item = psMetadataGetAndIncrement (iter)) != NULL) { 71 psMetadataAddItem (recipe, item, PS_LIST_TAIL, PS_META_REPLACE); 57 72 } 58 73 psFree (iter); 59 74 60 75 // set default recipe values here 61 METADATA_ADD_DEFAULT (config->recipe, Str, "FITMODE", "NONE", ""); 62 METADATA_ADD_DEFAULT (config->recipe, Str, "PHOTCODE", "NONE", ""); 63 METADATA_ADD_DEFAULT (config->recipe, Str, "BREAK_POINT", "NONE", ""); 64 METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_FORMAT", "CMP", ""); 65 METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_MODE", "SPLIT", ""); 76 // XXX this needs re-thinking... 77 METADATA_ADD_DEFAULT (recipe, Str, "FITMODE", "NONE", ""); 78 METADATA_ADD_DEFAULT (recipe, Str, "PHOTCODE", "NONE", ""); 79 METADATA_ADD_DEFAULT (recipe, Str, "BREAK_POINT", "NONE", ""); 66 80 67 81 // Chip selection: if we are using a single chip, select it for each FPA … … 88 102 } 89 103 } 90 91 // Construct cameras in preparation for reading92 input->fpa = pmFPAConstruct(config->camera);93 94 104 psTrace(__func__, 1, "Done with psphotParseCamera...\n"); 95 return input;105 return true; 96 106 } -
trunk/psphot/src/psphotReadout.c
r6522 r6571 2 2 3 3 // XXX 2006.02.07 : no leaks! 4 bool psphotReadout (pmReadout *readout, psMetadata *config) { 4 // XXX change 'config' to 'recipe' eventually 5 bool psphotReadout (pmConfig *config, pmFPAview *view) { 5 6 6 7 psArray *sources = NULL; 7 8 psArray *peaks = NULL; 8 psImage *skymodel = NULL;9 9 pmPSF *psf = NULL; 10 10 bool status; 11 11 12 // find the currently selected readout 13 psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); 14 pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); 15 pmReadout *readout = pmFPAviewThisReadout (view, input); 16 12 17 // generate a background model (median, smoothed image) 13 skymodel = psphotImageMedian (readout, config); 14 15 // optional output of background-subtracted image 16 psphot 18 psphotImageMedian (config, view); 17 19 18 20 // find the peaks in the image 19 peaks = psphotFindPeaks (readout, config);21 peaks = psphotFindPeaks (readout, recipe); 20 22 21 23 // construct sources and measure basic stats 22 24 // limit moments analysis by S/N? 23 sources = psphotSourceStats (readout, config, peaks);25 sources = psphotSourceStats (readout, recipe, peaks); 24 26 25 27 // classify sources based on moments, brightness 26 28 // faint sources not classified? 27 psphotRoughClass (sources, config);29 psphotRoughClass (sources, recipe); 28 30 29 31 // mark blended peaks PS_SOURCE_BLEND 30 psphotBasicDeblend (sources, config);32 psphotBasicDeblend (sources, recipe); 31 33 32 34 // use bright stellar objects to measure PSF 33 psf = psphotChoosePSF ( config, sources);35 psf = psphotChoosePSF (sources, recipe); 34 36 35 37 // linear PSF fit to peaks 36 psphotEnsemblePSF (readout, config, sources, psf, FALSE);38 psphotEnsemblePSF (readout, sources, recipe, psf, FALSE); 37 39 38 40 // non-linear PSF and EXT fit to brighter sources 39 psphotBlendFit (readout, config, sources, psf);41 psphotBlendFit (readout, sources, recipe, psf); 40 42 41 43 // replace fitted sources … … 46 48 47 49 // linear PSF fit to remaining peaks 48 psphotEnsemblePSF (readout, config, sources, psf, TRUE);50 psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); 49 51 50 52 // measure aperture photometry corrections 51 psphotApResid (readout, sources, config, psf);53 psphotApResid (readout, sources, recipe, psf); 52 54 53 55 // calculate source magnitudes 54 psphotMagnitudes ( config, sources, psf);56 psphotMagnitudes (sources, recipe, psf); 55 57 56 58 // update the header with stats results 59 // XXX need to do this conditionally? 57 60 psMetadata *header = pmReadoutGetHeader (readout); 58 61 psphotUpdateHeader (header, config); … … 60 63 // XXX make this an option? 61 64 // replace background in residual image 62 psphotSkyReplace ( readout, skymodel);65 psphotSkyReplace (config, view); 63 66 64 67 // need to do something with the sources, psf, and sky … … 67 70 // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.MEAN", PS_DATA_F32, "psphot sky mean", sky->sampleMean); 68 71 // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.SIGMA", PS_DATA_F32, "psphot sky stdev", sky->sampleStdev); 69 // XXX : what should I do with the skymodel image?70 72 71 73 // free up the local copies of the data … … 73 75 psFree (sources); 74 76 psFree (peaks); 75 psFree (skymodel);76 77 return true; 77 78 } -
trunk/psphot/src/psphotSkyReplace.c
r6481 r6571 1 1 # include "psphot.h" 2 2 3 // generate the median in NxN boxes, clipping heavily 4 // linear interpolation to generate full-scale model 5 bool psphotSkyReplace (pmReadout *readout, psImage *background) { 3 // in order to successfully replace the sky, we must define a corresponding file... 4 bool psphotSkyReplace (psMetadata *config, pmFPAview *view) { 6 5 6 // find the currently selected readout 7 psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); 8 pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); 9 pmReadout *readout = pmFPAviewThisReadout (view, input); 10 11 // select the corresponding images 7 12 psImage *image = readout->image; 8 13 psImage *mask = readout->mask; 14 15 // select background pixels, from output background file, or create 16 background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 0, 0, PS_TYPE_F32); 17 if (background == NULL) { 18 return false; 19 } 9 20 10 21 // replace the background model … … 17 28 } 18 29 30 psFree (background); 19 31 return true; 20 32 } -
trunk/psphot/src/testPSphotLoop.c
r6522 r6571 1 1 # include "psphot.h" 2 2 3 bool psphot Loop (psphotData *data, ppConfig *config) {3 bool psphotImageLoop (pmConfig *config) { 4 4 5 psRegion region = {0,0,0,0}; 6 pmFPA *fpa = data->input->fpa; 7 pmFPAiterator *fpi = pmFPAiteratorAlloc (fpa, region); 5 psMetadata *recipe = psMetadataLookupPts (NULL, config->recipes, PSPHOT_RECIPE); 6 if (!recipe) { 7 psErrorStackPrint(stderr, "Can't find recipe configuration!\n"); 8 exit(EXIT_FAILURE); 9 } 10 11 pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); 12 pmFPAview *view = pmFPAviewAlloc (0); 8 13 9 14 // files associated with the science image 10 pmFileDefine (fpi, config->camera, "PPIMAGE.INPUT"); 11 pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.HEAD"); 12 pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.DATA"); 13 pmFileDefine (fpi, config->camera, "PPIMAGE.BACKGND"); 14 pmFileDefine (fpi, config->camera, "PPIMAGE.BACKSUB"); 15 pmFileDefine (fpi, config->camera, "PPIMAGE.RESID"); 16 pmFileDefine (fpi, config->camera, "PPIMAGE.PSF"); 17 pmFileIOChecks (fpi); 15 pmFPAfileReadChecks (config->files, view); 18 16 19 while ((chip = pmChipNext ( fpi)) != NULL) {17 while ((chip = pmChipNext (view, input)) != NULL) { 20 18 21 19 psLogMsg ("psphot", 4, "Chip %d: %x %x\n", i, chip->exists, chip->process); 22 20 if (! chip->process) { continue; } 23 pmF ileIOChecks (fpi);21 pmFPAfileReadChecks (config->files, view); 24 22 25 while ((cell = pmCellNext ( fpi)) != NULL) {23 while ((cell = pmCellNext (view, input)) != NULL) { 26 24 27 25 psLogMsg ("psphot", 4, "Cell %d: %x %x\n", j, cell->exists, cell->process); 28 26 if (! cell->process) { continue; } 29 pmFileIOChecks (fpi); 27 pmFPAfileReadChecks (config->files, view); 28 29 // XXX optional mask and weight input image should be loaded here? 30 // this sets the weight map and basic mask applying CELL.BAD and CELL.SATURATION 31 pmCellSetWeights(cell); 32 33 // I have a valid mask, now mask in the analysis region of interest 34 pmCellSetMask (cell, recipe); 30 35 31 36 // process each of the readouts 32 while ((readout = pmReadoutNext ( fpi)) != NULL) {33 pmF ileIOChecks (fpi);37 while ((readout = pmReadoutNext (view, input)) != NULL) { 38 pmFPAfileReadChecks (config->files, view); 34 39 40 // run a single-model test if desired 41 // XXX move this to psphotReadout?? 42 // psphotModelTest (readout, recipe); 43 35 44 // run the actual photometry analysis 36 psphotReadout (readout, config); 45 psphotReadout (config, view); 46 47 pmFPAfileWriteChecks (config->files, view); 37 48 38 49 // write out the desired output dataset(s) 39 psphotOutput (readout, config);50 // psphotOutput (view, recipe); 40 51 } 52 pmFPAfileWriteChecks (config->files, view); 41 53 } 54 pmFPAfileWriteChecks (config->files, view); 42 55 } 43 pmF ileClose (fpi);56 pmFPAfileWriteChecks (config->files, view); 44 57 return true; 45 58 } 59 60 // I/O files related to psphot: 61 // PSPHOT.INPUT : input image file(s) 62 // PSPHOT.RESID : residual image 63 // PSPHOT.OUTPUT : output object tables (object) 64 65 // PSPHOT.BACKSUB : background subtracted image 66 // PSPHOT.BACKGND : background model (full-scale image?) 67 // PSPHOT.BACKMDL : background model (binned image?) 68 // PSPHOT.PSF : sample PSF images 69 70 71 /** 72 73 filename | pmFPAfile | pmFPA | pmFPAview 74 chip00.f | PSPHOT.IN | input | view 75 chip01.f | PSPHOT.IN | input | view 76 chip02.f | PSPHOT.IN | input | view 77 chip03.f | PSPHOT.IN | input | view 78 79 out00.f | PSPHOT.OUT | input | view 80 out01.f | PSPHOT.OUT | input | view 81 out02.f | PSPHOT.OUT | input | view 82 out03.f | PSPHOT.OUT | input | view 83 84 obj00.f | PSPHOT.OBJ | input | view 85 obj01.f | PSPHOT.OBJ | input | view 86 obj02.f | PSPHOT.OBJ | input | view 87 obj03.f | PSPHOT.OBJ | input | view 88 89 bin00.f | PSPHOT.BIN | binned | view 90 bin01.f | PSPHOT.BIN | binned | view 91 bin02.f | PSPHOT.BIN | binned | view 92 bin03.f | PSPHOT.BIN | binned | view 93 94 mosaic.f | PSPHOT.MOS | mosaic | view 95 96 **/
Note:
See TracChangeset
for help on using the changeset viewer.
