Index: anches/eam_rel9_b1/psModules/src/objects/psLibUtils.c =================================================================== --- /branches/eam_rel9_b1/psModules/src/objects/psLibUtils.c (revision 5854) +++ (revision ) @@ -1,509 +1,0 @@ -# include // for strncasecmp -# include -# include "psLibUtils.h" - -// XXX EAM : this is NOT included in the C99 headers ?? -FILE *fdopen(int fildes, const char *mode); - -// XXX EAM : these utility functions should be added back into PSLib - -static psHash *timers = NULL; - -/* XXX: remove -bool psTimerClear (char *name) { - - bool status; - - if (name == NULL) return false; - - status = psHashRemove (timers, name); - return (status); -} -*/ - -void psTimerFree () -{ - - psFree (timers); - p_psTimeFinalize(); - return; -} - -// start/restart a named timer -bool psTimerStart (char *name) -{ - - psTime *start; - - if (timers == NULL) - timers = psHashAlloc (16); - - start = psTimeGetNow (PS_TIME_UTC); - psHashAdd (timers, name, start); - psFree (start); - return (TRUE); -} - -// get current elapsed time on named timer -psF64 psTimerMark (char *name) -{ - - psTime *start; - psTime *mark; - psF64 delta; - - if (timers == NULL) - return (0); - - start = psHashLookup (timers, name); - if (start == NULL) - return (0); - - mark = psTimeGetNow (PS_TIME_UTC); - delta = psTimeDelta (mark, start); - psFree (mark); - // psFree (start); -- XXX is psHashLookup not psMemCopying? - - return (delta); -} - -/* XXX: remove -// find the location of the specified argument -int psArgumentGet (int argc, char **argv, char *arg) { - - int i; - - for (i = 0; i < argc; i++) { - if (!strcmp(argv[i], arg)) - return (i); - } - - return ((int)NULL); -} -*/ - -/* XXX: remove -// remove the specified argument (by location) -int psArgumentRemove (int N, int *argc, char **argv) { - - int i; - - if ((N != (int)NULL) && (N != 0)) { - (*argc)--; - for (i = N; i < *argc; i++) { - argv[i] = argv[i+1]; - } - } - - return (N); -} -*/ - -// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) -// 2 = default, -v = 3, -vv = 4, -vvv = 5 -psS32 psLogArguments (int *argc, char **argv) -{ - - int N, level; - - // default log level is 2 - level = 2; - - // set in order, so that -vvv overrides -vv overrides -v - if ((N = psArgumentGet (*argc, argv, "-v"))) { - psArgumentRemove (N, argc, argv); - level = 3; - } - if ((N = psArgumentGet (*argc, argv, "-vv"))) { - psArgumentRemove (N, argc, argv); - level = 4; - } - if ((N = psArgumentGet (*argc, argv, "-vvv"))) { - psArgumentRemove (N, argc, argv); - level = 5; - } - - if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { - if (*argc < N + 2) { - psAbort ("psLogArguments", "USAGE: -logfmt (format)"); - } - psArgumentRemove (N, argc, argv); - psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid - psArgumentRemove (N, argc, argv); - } - - // set the level, return the level - psLogSetLevel (level); - return (level); -} - -// set trace levels by facility -psS32 psTraceArguments (int *argc, char **argv) -{ - - int N; - - // argument format is: -trace (facil) (level) - while ((N = psArgumentGet (*argc, argv, "-trace"))) { - if (*argc < N + 3) { - psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); - } - psArgumentRemove (N, argc, argv); - psTraceSetLevel (argv[N], atoi(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { - psTracePrintLevels (); - exit (2); - } - return (TRUE); -} - -# if (0) - // alternate implementation of this function from pmObjects.c - // now defined in psLib SDRS as psImageRow - psVector *psGetRowVectorFromImage(psImage *image, - psU32 row) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32)); - return(tmpVector); -} -# endif - -// XXX EAM : this is now in psLib -void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma) -{ - - int Nx, Ny, Npixel, Nrange; - float factor, g, s; - psVector *temp; - - // relevant terms - Nrange = sigma*Nsigma + 0.5; - Npixel = 2*Nrange + 1; - factor = -0.5/(sigma*sigma); - - Nx = image->numCols; - Ny = image->numRows; - - // generate gaussian - psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32); - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] = exp (factor*i*i); - } - psF32 *gauss = &gaussnorm->data.F32[Nrange]; - - // smooth in X direction - temp = psVectorAlloc (Nx, PS_TYPE_F32); - for (int j = 0; j < Ny; j++) { - psF32 *vi = image->data.F32[j]; - psF32 *vo = temp->data.F32; - for (int i = 0; i < Nx; i++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (i+n < 0) - continue; - if (i+n >= Nx) - continue; - s += gauss[n]*vi[i+n]; - g += gauss[n]; - } - vo[i] = s / g; - } - memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32)); - } - psFree (temp); - - // smooth in Y direction - temp = psVectorAlloc (image->numRows, PS_TYPE_F32); - for (int i = 0; i < Nx; i++) { - psF32 *vo = temp->data.F32; - psF32 **vi = image->data.F32; - for (int j = 0; j < Ny; j++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (j+n < 0) - continue; - if (j+n >= Ny) - continue; - s += gauss[n]*vi[j+n][i]; - g += gauss[n]; - } - vo[j] = s / g; - } - // replace temp in image - for (int j = 0; j < Ny; j++) { - vi[j][i] = vo[j]; - } - } - psFree (temp); - psFree (gaussnorm); -} - -bool psImageInit (psImage *image,...) -{ - - va_list argp; - psU8 vU8; - psF32 vF32; - psF64 vF64; - - if (image == NULL) - return (false); - - va_start (argp, image); - - switch (image->type.type) { - case PS_TYPE_U8: - vU8 = va_arg (argp, psU32); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.U8[iy][ix] = vU8; - } - } - break; - - case PS_TYPE_F32: - vF32 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F32[iy][ix] = vF32; - } - } - return (true); - - case PS_TYPE_F64: - vF64 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F64[iy][ix] = vF64; - } - } - return (true); - - default: - psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type); - return (false); - } - return (false); -} - -/* XXX: remove -// count number of pixels with given mask value -int psImageCountPixelMask (psImage *mask, psU8 value) -{ - int Npixels = 0; - - for (int i = 0; i < mask->numRows; i++) { - for (int j = 0; j < mask->numCols; j++) { - if (mask->data.U8[i][j] & value) { - Npixels ++; - } - } - } - return (Npixels); -} -*/ - -// define a square region centered on the given coordinate -// XXX EAM : this is now in psLib -# if (0) - psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) -{ - psRegion *region; - region = psRegionAlloc (x - radius, x + radius + 1, - y - radius, y + radius + 1); - return (region); -} -# endif - -// set actual region based on image parameters: -// compensate for negative upper limits -// XXX this is inconsistent: the coordindates should always be in the parent -// frame, which means the negative values should subtract from Nx,Ny of -// the parent, not the child. but, we don't carry the dimensions of the -// parent in the psImage container. for now, use the child Nx,Ny -// force range to be on this subimage -// XXX EAM : this needs to be changed to use psRegion rather than psRegion* -// XXX EAM : this is now in psLib -# if (0) - psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) -{ - - // x0,y0, x1,y1 are in *parent* units - - if (out == NULL) { - out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1); - } else { - *out = *in; - } - // XXX these are probably wrong (see above) - if (out->x1 <= 0) { - out->x1 = image->col0 + image->numCols + out->x1; - } - if (out->y1 <= 0) { - out->y1 = image->row0 + image->numRows + out->y1; - } - - // force the lower-limits to be on the child - out->x0 = PS_MAX(image->col0, out->x0); - out->y0 = PS_MAX(image->row0, out->y0); - - // force the upper-limits to be on the child - out->x1 = PS_MIN(image->col0 + image->numCols, out->x1); - out->y1 = PS_MIN(image->row0 + image->numRows, out->y1); - return (out); -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) -{ - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) - continue; - if (ix + image->col0 >= region->x1) - continue; - if (iy + image->row0 < region->y0) - continue; - if (iy + image->row0 >= region->y1) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area not contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) -{ - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) - goto maskit; - if (ix + image->col0 >= region->x1) - goto maskit; - if (iy + image->row0 < region->y0) - goto maskit; - if (iy + image->row0 >= region->y1) - goto maskit; - continue; -maskit: - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) -{ - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 > R2) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) -{ - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 < R2) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -/* XXX: remove -psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { - - int nBin = (upper - lower) / delta; - - psVector *out = psVectorAlloc (nBin, type); - - for (int i = 0; i < nBin; i++) { - out->data.F64[i] = lower + i * delta; - } - - return (out); -} -*/ - -// XXX EAM a utility function -bool p_psVectorPrintRow (int fd, psVector *a, char *name) -{ - - FILE *f; - f = fdopen(fd, "a+"); - fprintf (f, "vector: %s\n", name); - - for (int i = 0; i < a[0].n; i++) { - fprintf (f, "%f ", p_psVectorGetElementF64(a, i)); - } - fprintf (f, "\n"); - fclose(f); - return (true); -} -// XXX EAM is the use of fdopen here a good way to do this? Index: anches/eam_rel9_b1/psModules/src/objects/psLibUtils.h =================================================================== --- /branches/eam_rel9_b1/psModules/src/objects/psLibUtils.h (revision 5854) +++ (revision ) @@ -1,59 +1,0 @@ - -# ifndef PS_LIB_UTILS -# define PS_LIB_UTILS - -// XXX EAM : psEllipse needs to be be included in psLib -# include "psEllipse.h" - -// structure to carry a dynamic string -typedef struct -{ - int NLINE; - int Nline; - char *line; -} -psLine; - -# define psMemCopy(A)(psMemIncrRefCounter((A))) - -// XXX EAM : my version using varience instead of stdev -bool psMinimizeGaussNewtonDelta_EAM (psVector *delta, - const psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func); - -// minimize : using varience vs sigma and parameter limits -psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, psVector *Beta, psVector *Params, const psImage *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda); -psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func); -psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda); - -// psLib extra utilities -bool psTimerStart (char *name); // added to SDRS -//XXX: I removed this: bool psTimerClear (char *name); // added to SDRS -psF64 psTimerMark (char *name); // added to SDRS -void psTimerFree (); // added to SDRS (as psTimerStop) -psS32 psLogArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -psS32 psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -//XXX: remove: int psArgumentGet (int argc, char **argv, char *arg); // added to SDRS -//XXX: remove: int psArgumentRemove (int N, int *argc, char **argv); // added to SDRS -//XXX: remove: psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS -// psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow) -//XXX: remove: int psImageCountPixelMask (psImage *mask, psU8 value); // added to SDRS - -// basic image functions -bool psImageInit (psImage *image,...); // added to SDRS (v.15) -void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma); // added to SDRS (v.15) - -// psLine functions -- keep out for now? -psLine *psLineAlloc (int Nline); -bool psLineInit (psLine *line); -bool psLineAdd (psLine *line, char *format, ...); - -// not included in the .h file -- these are fairly lame, keep out? -bool p_psVectorPrint (int fd, psVector *a, char *name); -bool p_psVectorPrintRow (int fd, psVector *a, char *name); - -# endif Index: anches/eam_rel9_b1/psModules/src/objects/psModulesUtils.c =================================================================== --- /branches/eam_rel9_b1/psModules/src/objects/psModulesUtils.c (revision 5854) +++ (revision ) @@ -1,95 +1,0 @@ -# include "psModulesUtils.h" - -// XXXX: Must figure out the right names for these defines, then replace! -#define PS_META_STR 0 -#define PS_META_F32 0 - - - -// extract config informatin from config data or from image header, if specified -// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) -{ - - char *source; - char *keyword; - psF32 value; - psMetadataItem *item; - - // look for the entry in the config collection - item = psMetadataLookup (config, name); - if (item == NULL) { - psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name); - value = psMetadataLookupF32 (status, header, name); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in header", name); - } - *status = true; - return (value); - } - - // I'm either expecting a string, with the name "HD:keyword"... - if (item->type == PS_META_STR) { - source = item->data.V; - if (!strncasecmp (source, "HD:", 3)) { - keyword = &source[3]; - value = psMetadataLookupF32 (status, header, keyword); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in config", name); - } - *status = true; - // psFree (item); - return (value); - } - } - - // ... or a value (F32?) - if (item->type == PS_META_F32) { - value = item->data.F32; - // psFree (item); - return (value); - } - - *status = false; - psError(PS_ERR_UNKNOWN, true, "invalid item"); - return (0); -} - -bool pmModelFitStatus (pmModel *model) -{ - - bool status; - - pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); - status = statusFunc (model); - - return (status); -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) -{ - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { - if (mask->data.U8[iy][ix]) - continue; - obsSum += image->data.F32[iy][ix] - sky; - } - } - if (obsSum <= 0) - return false; - if (fitSum <= 0) - return false; - - *fitMag = -2.5*log10(fitSum); - *obsMag = -2.5*log10(obsSum); - return (true); -} - Index: anches/eam_rel9_b1/psModules/src/objects/psModulesUtils.h =================================================================== --- /branches/eam_rel9_b1/psModules/src/objects/psModulesUtils.h (revision 5854) +++ (revision ) @@ -1,20 +1,0 @@ - -#ifndef PS_MODULE_UTILS -#define PS_MODULE_UTILS - -#include // for strcasecmp -#include -#include "psLibUtils.h" -#include "pmObjects.h" -#include "pmModelGroup.h" - -// psModule extra utilities -// bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); -bool pmModelFitStatus (pmModel *model); -int pmSourceDophotType (pmSource *source); -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask); - -// XXX: unify with paul's image/header/metadata functions -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); - -#endif