Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 4953) +++ /trunk/psphot/Makefile (revision 4954) @@ -14,6 +14,7 @@ INCS = $(IPSLIB) -LIBS = -lpsmodule $(LPSLIB) -CFLAGS = $(INCS) -std=c99 +# LIBS = -lpsmodule $(LPSLIB) +LIBS = $(LPSLIB) +CFLAGS = $(INCS) -std=c99 -Wall -Werror LFLAGS = $(LIBS) @@ -30,4 +31,5 @@ $(SRC)/psphotMarkPSF.$(ARCH).o \ $(SRC)/psphotSubtractPSF.$(ARCH).o \ +$(SRC)/psphotSortBySN.$(ARCH).o \ $(SRC)/pspsf.$(ARCH).o \ $(SRC)/psLibUtils.$(ARCH).o \ @@ -39,5 +41,7 @@ $(SRC)/psMinimize.$(ARCH).o \ $(SRC)/psImageData.$(ARCH).o \ -$(SRC)/psphotSortBySN.$(ARCH).o +$(SRC)/psEllipse.$(ARCH).o \ +$(SRC)/pmModelInit.$(ARCH).o \ +$(SRC)/pmObjects_EAM.$(ARCH).o FITSOURCE = \ Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 4953) +++ /trunk/psphot/doc/notes.txt (revision 4954) @@ -1,4 +1,64 @@ Notes on psphot + +2005.09.06 + + I have built a working version of PSPhot using my own copy of the + pmObjects.[ch] files. I need to minimize the difference between the + v0.7.0, v0.5.0, and my own version of pmObjects.c + + changes v0.5.0 / v0.7.0 / EAM + - changed all ps* to pm* + - fixed naming for local static functions (eg modelFree not p_psModelFree) + - dropped all hard-coded model names (vs v0.7.0) + - using current psMemSetDellocator name + - added noise image to pmSource + - converted old asserts to new asserts (PS_.._CHECK..) + - add return from error in pmFindImagePeaks + - allow isItInThisRegion to have NULL region (true) + - fixed subImages to correspond to region definition (end is + exclusive) in LocalSky and SetPixelCircle + - dropped old, redundant code in pmSourceLocalSky + - fixed pmSourceMoments use of mask (0 is valid, not 1) + - added pmSourceMoments validity tests (numPixels, Sum, centroid shift) + - added correction for Sxy (XY/Sum - x*y) + - added pmSourcePSFClump function (and pmPSFClump structure) + - allow missing sources in pmSourcePSFClump (not all peaks yield + sources) + - moved pmSourceRoughClass metadata lookups out of inner loop + - added saturated pixel test using pmImageCountPixelMask + - various redefinitions of initial classes + - dropped old, redundant code in pmSourceSetPixelsCircle + - mask image set to type psU8 (was psF32) + - abstract object model functions for pmSourceModelGuess + - converted evalModel to a public function psModelEval, modified params + - pmSourceFitModel operates on a specified model + - pmSourceFitModel had option to apply/ignore PSF information + - pmSourceFitModel calculates yErr + - pmSourceFitModel calculates the covariance matrix + - added paramMask to pmSourceFitModel + - pmSourceFitModel_EAM uses covar to carry in beta and param + limits... + - pmSourceFitModel requires solution to stay within image + - pmSourceFitModel calculates and saves chisq, nIter, nDOF + - pmSourceFitModel uses GaussNewtonDelta for frozen params + - added mask to pmSourceAdd, pmSourceSub + - model abstractions in pmSourceAdd, pmSourceSub + - fixed 'center' option in pmSourceAdd/Sub + + - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var)) + - pmSourceContour is using the model flux, not the image flux + - does psImageAlloc zero the image? if so, drop init routine + from pmSourcePSFClump + - pmSourceRoughClass uses RDNOISE, GAIN to calculate SN: should use + the source->noise image? part of pmSourceMoments? + - pmSourceMoments_EAM uses macro for checkRegion + - pmObjects_EAM.c uses some hardcoded mask values (LocalSky) + - why does LocalSky_EAM not need subImageMask->col0 = subImage->col0?? + - psGetRowVectorFromImage should be in psLib (in psLibUtils) + - pmSourceSetPixelsCircle and pmSourceLocalSky probably could use + the psImageMaskRegion, etc functions. + - pmSourceSetPixelsCircle uses a hard-coded mask value + - 2005.09.05 Index: /trunk/psphot/doc/psphot.txt =================================================================== --- /trunk/psphot/doc/psphot.txt (revision 4953) +++ /trunk/psphot/doc/psphot.txt (revision 4954) @@ -1,5 +1,38 @@ -psphot is a stand-alone program which implements the Pan-STARRS IPP -object detection and analysis functions for using with individual -images. +Defined APIs: +pmMoments *pmMomentsAlloc(); +pmModel *pmModelAlloc(); +pmSource *pmSourceAlloc(); +psVector *pmFindVectorPeaks() +psArray *pmFindImagePeaks() +psList *pmCullPeaks() +pmSource *pmSourceLocalSky() +bool pmSourceMoments() +pmPSFClump pmSourcePSFClump() +bool pmSourceRoughClass() +bool pmSourceSetPixelsCircle() +pmModel *pmSourceModelGuess() +psArray *pmSourceContour() +bool pmSourceFitModel() +bool pmSourceAddModel() +bool pmSourceSubModel() +int pmModelParameterCount (); +char *pmModelGetType (); +pmModelType pmModelSetType (); +pmModelFunc pmModelFunc_GetFunction (); +pmModelFlux pmModelFlux_GetFunction (); +pmModelRadius pmModelRadius_GetFunction (); +pmModelLimits pmModelLimits_GetFunction (); +pmModelGuessFunc pmModelGuessFunc_GetFunction (); +pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (); +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (); + +model function abstractions: + typedef psMinimizeLMChi2Func pmModelFunc; + typedef psF64 (*pmModelFlux)(const psVector *params); + typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); + typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); + typedef bool (*pmModelFitStatusFunc)(pmModel *model); Index: unk/psphot/src/fitsource.c =================================================================== --- /trunk/psphot/src/fitsource.c (revision 4953) +++ (revision ) @@ -1,21 +1,0 @@ -# include "psphot.h" - -int main (int argc, char **argv) { - - psMetadata *config = NULL; // user-provided configuration information - psImageData *imdata = NULL; - float x, y; - - config = fs_args (&argc, argv); - - // load input data (image and config) - // create or load mask and noise images - imdata = setup (config); - - x = atof (argv[4]); - y = atof (argv[5]); - - onesource (imdata, config, argv[3], x, y); - - exit (0); -} Index: unk/psphot/src/fs_args.c =================================================================== --- /trunk/psphot/src/fs_args.c (revision 4953) +++ (revision ) @@ -1,49 +1,0 @@ -# include "psphot.h" - -psMetadata *fs_args (int *argc, char **argv) { - - int N, Nfail; - - // basic pslib options - psLogSetFormat ("M"); - psLogArguments (argc, argv); - psTraceArguments (argc, argv); - - // identify options in args list - these go on config - - char *mask = NULL; - if ((N = get_argument (*argc, argv, "--mask"))) { - remove_argument (N, argc, argv); - mask = psStringCopy (argv[N]); - remove_argument (N, argc, argv); - } - - char *noise = NULL; - if ((N = get_argument (*argc, argv, "--noise"))) { - remove_argument (N, argc, argv); - noise = psStringCopy (argv[N]); - remove_argument (N, argc, argv); - } - - if (*argc != 6) fs_usage (); - - psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[2], FALSE); - - psMetadataAddStr (config, PS_LIST_HEAD, "INPUT", "", argv[1]); - if (mask != NULL) { - psMetadataAddStr (config, PS_LIST_HEAD, "MASK", "", mask); - } - if (noise != NULL) { - psMetadataAddStr (config, PS_LIST_HEAD, "NOISE", "", noise); - } - return (config); -} - -int fs_usage () { - - fprintf (stderr, "USAGE: psphot (image.fits) (config) (model) (x) (y)\n"); - fprintf (stderr, "options: \n"); - fprintf (stderr, " --mask (filename)\n"); - fprintf (stderr, " --noise (filename)\n"); - exit (2); -} Index: /trunk/psphot/src/models/pmModel_GAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_GAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_GAUSS.c (revision 4954) @@ -0,0 +1,152 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF32 pmModelFunc_GAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 r = exp(-z); + psF32 q = params->data.F32[1]*r; + psF32 f = q + params->data.F32[0]; + + if (deriv != NULL) { + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +psF64 pmModelFlux_GAUSS(const psVector *params) +{ + psF64 A1 = PS_SQR(params->data.F32[4]); + psF64 A2 = PS_SQR(params->data.F32[5]); + psF64 A3 = PS_SQR(params->data.F32[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 pmModelRadius_GAUSS (const psVector *params, psF64 flux) +{ + + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1]) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +// define the parameter limits +bool pmModelLimits_GAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + + return (TRUE); +} + +// make an initial guess for parameters +bool pmModelGuess_GAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = moments->x; + params[3] = moments->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + return(true); +} + +// construct the PSF model from the FLT model and the psf +bool pmModelFromPSF_GAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +// check the status of the fitted model +bool pmModelFitStatus_GAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /trunk/psphot/src/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_PGAUSS.c (revision 4954) @@ -0,0 +1,167 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF32 pmModelFunc_PGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + psF32 t = 1 + z + z*z/2.0; + psF32 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = PAR[1]*r + psF32 q = PAR[1]*r*r*t; + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +bool pmModelLimits_PGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (7, PS_TYPE_F32); + *params_min = psVectorAlloc (7, PS_TYPE_F32); + *params_max = psVectorAlloc (7, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + + return (TRUE); +} + +psF64 pmModelFlux_PGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + z + z*z/2 + z*z*z/6); + norm += f; + } + norm *= 0.01; + + psF64 Flux = params->data.F32[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_PGAUSS (const psVector *params, psF64 flux) +{ + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1]) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = moments->x; + params[3] = moments->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + return(true); +} + +bool pmModelFromPSF_PGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 7; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_PGAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /trunk/psphot/src/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_QGAUSS.c (revision 4954) @@ -0,0 +1,192 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^M + z^N) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = + params->data.F32[8] = +*****************************************************************************/ + +psF32 pmModelFunc_QGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + psF32 r = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25)); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*z; + } + return(f); +} + +bool pmModelLimits_QGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (8, PS_TYPE_F32); + *params_min = psVectorAlloc (8, PS_TYPE_F32); + *params_max = psVectorAlloc (8, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + beta_lim[0][0].data.F32[7] = 0.5; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + params_min[0][0].data.F32[7] = 0.1; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +5.0; + params_max[0][0].data.F32[7] = 10.0; + + return (TRUE); +} + +bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source) { + + pmMoments *moments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + params[0] = moments->Sky; + params[1] = moments->Peak - moments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.2 / moments->Sx; + params[5] = 1.2 / moments->Sy; + params[6] = 0.0; + params[7] = 1.0; + return(true); +} + +psF64 pmModelFlux_QGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_QGAUSS (const psVector *params, psF64 flux) +{ + psF64 z, f; + psF32 *PAR = params->data.F32; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // if Sx == Sy, sigma = Sx == Sy + psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + psF64 dz = 1.0 / (2.0 * sigma*sigma); + psF64 limit = flux / PAR[1]; + + // we can do this much better with intelligent choices here + for (z = 0.0; z < 20.0; z += dz) { + f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + if (f < limit) break; + } + psF64 radius = sigma * sqrt (2.0 * z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelFromPSF_QGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_QGAUSS (pmModel *model) { + + psF32 dP; + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + + if (status) return true; + return false; +} Index: /trunk/psphot/src/models/pmModel_RGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_RGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_RGAUSS.c (revision 4954) @@ -0,0 +1,151 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z + z^Npow) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow +*****************************************************************************/ + +psF64 psModelFunc_RGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + // psF32 FACT = 1 + 5*exp(-5*PAR[7]); + + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z + z*p); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(1 + PAR[7]*p); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -q*px*X*2; + deriv->data.F32[5] = -q*py*Y*2; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -5.0*t*log(z)*p*z; + + // deriv->data.F32[4] = -1.8*q*px*X*2; + // deriv->data.F32[5] = -1.8*q*py*Y*2; + // deriv->data.F32[6] = -1.5*q*X*Y; + // deriv->data.F32[7] = -5.0*t*log(z)*p*z; + } + return(f); +} + +psF64 psModelFlux_RGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + z + pow(z, PAR[7])); + norm += f; + } + norm *= 0.01; + + psF64 Flux = params->data.F32[1] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 psModelRadius_RGAUSS (const psVector *params, psF64 flux) +{ + psF64 z, f, p; + psF32 *PAR = params->data.F32; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // if Sx == Sy, sigma = Sx == Sy + psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0); + psF64 dz = 1.0 / (2.0 * sigma*sigma); + psF64 limit = flux / PAR[1]; + + // we can do this much better with intelligent choices here + for (z = 0.0; z < 20.0; z += dz) { + p = pow(z, PAR[7]); + f = 1.0 / (1 + z + p); + if (f < limit) break; + } + psF64 radius = sigma * sqrt (2.0 * z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool psModelGuess_RGAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(source->moments->Sx); + moments.y2 = PS_SQR(source->moments->Sy); + moments.xy = source->moments->Sxy; + + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = 1.0 / shape.sx; + params->data.F32[5] = 1.0 / shape.sy; + params->data.F32[6] = shape.sxy; + params->data.F32[7] = 2.0; + return(true); +} + +bool psModelFromPSF_RGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /trunk/psphot/src/models/pmModel_SGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_SGAUSS.c (revision 4954) @@ -0,0 +1,342 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^Npow + St*z^Ntide) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow + params->data.F32[8] = St +*****************************************************************************/ + +# define SQ(A)((A)*(A)) +psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image); +psF64 p_psImageGetElementF64(psImage *a, int i, int j); + +psF32 pmModelFunc_SGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8); + // note that if z -> 0, dPAR[7] -> -inf + // also z^(PAR[7]-1) -> Inf + + psF32 pr = z*PAR[8]; + psF32 pr3 = pr*pr*pr; + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z*p + pr*pr3); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -2*t*log(z)*z*p; + deriv->data.F32[8] = -2*t*4*z*pr3; + } + return(f); +} + +bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) { + + *beta_lim = psVectorAlloc (9, PS_TYPE_F32); + *params_min = psVectorAlloc (9, PS_TYPE_F32); + *params_max = psVectorAlloc (9, PS_TYPE_F32); + + beta_lim[0][0].data.F32[0] = 1000; + beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[2] = 5; + beta_lim[0][0].data.F32[3] = 5; + beta_lim[0][0].data.F32[4] = 0.5; + beta_lim[0][0].data.F32[5] = 0.5; + beta_lim[0][0].data.F32[6] = 0.5; + beta_lim[0][0].data.F32[7] = 0.5; + beta_lim[0][0].data.F32[8] = 0.05; + + params_min[0][0].data.F32[0] = -1000; + params_min[0][0].data.F32[1] = 0; + params_min[0][0].data.F32[2] = -100; + params_min[0][0].data.F32[3] = -100; + params_min[0][0].data.F32[4] = 0.01; + params_min[0][0].data.F32[5] = 0.01; + params_min[0][0].data.F32[6] = -5.0; + params_min[0][0].data.F32[7] = 0.5; + params_min[0][0].data.F32[8] = 0.001; + + params_max[0][0].data.F32[0] = 1e5; + params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! + params_max[0][0].data.F32[4] = 2.0; + params_max[0][0].data.F32[5] = 2.0; + params_max[0][0].data.F32[6] = +3.0; + params_max[0][0].data.F32[7] = 5.0; + params_max[0][0].data.F32[8] = 0.5; + + return (TRUE); +} + +// measure the flux for the elliptical contour +psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image) { + + double t, dt, ct, st, xo, yo, value; + int N, Nt, x, y; + + // choose dt to uniformly divide contour, with ~1 pix spacing at most + dt = asin (1 / axes.minor); + Nt = (int)(2*M_PI / dt) + 1; + dt = 2*M_PI / Nt; + + ct = cos(axes.theta); + st = sin(axes.theta); + xo = xc - image->col0; + yo = yc - image->row0; + + psVector *contour = psVectorAlloc (Nt, PS_TYPE_F32); + for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) { + x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo; + y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo; + value = p_psImageGetElementF64(image, x, y); + if (isfinite(value)) { + contour->data.F32[N] = value; + N++; + } + } + contour->n = N; + // accept every pixel: double counting is not so problematic here... + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, contour, NULL, NULL, 0); + value = stats->sampleMedian; + + psFree (stats); + psFree (contour); + + return (value); +} + +bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source) { + + pmMoments *sMoments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(sMoments->Sx); + moments.y2 = PS_SQR(sMoments->Sy); + moments.xy = sMoments->Sxy; + + // solve the math to go from Moments To Shape + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params[0] = sMoments->Sky; + params[1] = sMoments->Peak - sMoments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + +# if (0) + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + axes.major *= 2.0; + axes.minor *= 2.0; + f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + + if (f1 > f2) { + params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); + } else { + params[7] = 0.5; + } +# endif + + params[7] = 1.8; + params[8] = 0.1; + + + return(true); +} + +// XXX EAM : test version using flux contours to guess slope +bool pmModelGuess_SGAUSS_HARD (pmModel *model, pmSource *source) { + + pmMoments *sMoments = source->moments; + pmPeak *peak = source->peak; + psF32 *params = model->params->data.F32; + float f1, f2; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(sMoments->Sx); + moments.y2 = PS_SQR(sMoments->Sy); + moments.xy = sMoments->Sxy; + + // solve the math to go from Moments To Shape + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params[0] = sMoments->Sky; + params[1] = sMoments->Peak - sMoments->Sky; + params[2] = peak->x; + params[3] = peak->y; + params[4] = 1.0 / shape.sx; + params[5] = 1.0 / shape.sy; + params[6] = shape.sxy; + + f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + axes.major *= 2.0; + axes.minor *= 2.0; + f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); + + if (f1 > f2) { + params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); + } else { + params[7] = 0.5; + } + params[8] = 0.1; + + return(true); +} + +psF64 pmModelFlux_SGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + for (z = 0.005; z < 50; z += 0.01) { + psF32 pr = PAR[8]*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 pmModelRadius_SGAUSS (const psVector *params, psF64 flux) +{ + psF64 r, z, pr, f; + psF32 *PAR = params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // convert Sx,Sy,Sxy to major/minor axes + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = SQ(r); + pr = PAR[8]*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + if (f < limit) break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool pmModelFromPSF_SGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 9; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} + +bool pmModelFitStatus_SGAUSS (pmModel *model) { + + psF32 dP; + bool status; + EllipseAxes axes; + EllipseShape shape; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + + dP = 0; + dP += PS_SQR(dPAR[4] / PAR[4]); + dP += PS_SQR(dPAR[5] / PAR[5]); + dP += PS_SQR(dPAR[7] / PAR[7]); + dP = sqrt (dP); + + status = true; + status &= (dP < 0.5); + status &= (PAR[1] > 0); + status &= ((dPAR[1]/PAR[1]) < 0.5); + status &= (fabs(PAR[8]) < 0.5); + status &= (dPAR[8] < 0.1); + status &= (axes.major > 1.41); + status &= (axes.minor > 1.41); + status &= ((axes.major / axes.minor) < 5.0); + status &= (PAR[7] > 0.5); + + if (status) return true; + return false; +} Index: /trunk/psphot/src/models/pmModel_TGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_TGAUSS.c (revision 4954) @@ -0,0 +1,109 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^M + z^N) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = +*****************************************************************************/ + +psF64 psModelFunc_TGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + psF32 p = pow(z, 1.2); + psF32 r = 1.0 / (1 + z + PAR[7]*z*p); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(1 + PAR[7]*2.2*p); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*z*p; + } + return(f); +} + +psF64 psModelFlux_TGAUSS(const psVector *params) +{ + psF64 A1 = 1 / PS_SQR(params->data.F32[4]); + psF64 A2 = 1 / PS_SQR(params->data.F32[5]); + psF64 A3 = params->data.F32[6]; + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 psModelRadius_TGAUSS (const psVector *params, psF64 flux) +{ + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1]) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool psModelGuess_TGAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = 1.0/source->moments->Sx; + params->data.F32[5] = 1.0/source->moments->Sy; + // params->data.F32[6] = source->moments->Sxy; + params->data.F32[6] = 0.0; + params->data.F32[7] = 5.0; + return(true); +} + +bool psModelFromPSF_TGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /trunk/psphot/src/models/pmModel_WAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_WAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_WAUSS.c (revision 4954) @@ -0,0 +1,100 @@ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ + +psF64 psModelFunc_WAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[2]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); + psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + if (deriv != NULL) { + // note difference from gaussian: q = params->data.F32[1]*r + psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; + deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; + // The values of 100 dampen the swing of params->data.F32[7,8] */ + } + return(f); +} + +// this is probably wrong since it uses the gauss integral 2 pi sigma^2 +psF64 psModelFlux_WAUSS(const psVector *params) +{ + psF64 A1 = 1 / PS_SQR(params->data.F32[4]); + psF64 A2 = 1 / PS_SQR(params->data.F32[5]); + psF64 A3 = params->data.F32[6]; + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); + // Area is equivalent to 2 pi sigma^2 + + psF64 Flux = params->data.F32[1] * Area; + + return(Flux); +} + +// return the radius which yields the requested flux +psF64 psModelRadius_WAUSS (const psVector *params, psF64 flux) +{ + if (flux <= 0) return (1.0); + if (params->data.F32[1] <= 0) return (1.0); + if (flux >= params->data.F32[1] <= 0) return (1.0); + + psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); + psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); + return (radius); +} + +bool psModelGuess_WAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + // XXX: What are these? + // params->data.F32[7] = B2; + // params->data.F32[8] = B3; + return(true); +} + +bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 9; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: /trunk/psphot/src/models/pmModel_ZGAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_ZGAUSS.c (revision 4954) +++ /trunk/psphot/src/models/pmModel_ZGAUSS.c (revision 4954) @@ -0,0 +1,158 @@ + +/****************************************************************************** + one component, two slopes: + 1 / (1 + z^Npow + PAR8*z^4) + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = 1 / SigmaX; + params->data.F32[5] = 1 / SigmaY; + params->data.F32[6] = Sxy; + params->data.F32[7] = Npow +*****************************************************************************/ + +# define SQ(A)((A)*(A)) +# define PAR8 0.1 + +psF64 psModelFunc_ZGAUSS(psVector *deriv, + const psVector *params, + const psVector *x) +{ + psF32 *PAR = params->data.F32; + + psF32 X = x->data.F32[0] - PAR[2]; + psF32 Y = x->data.F32[1] - PAR[3]; + psF32 px = PAR[4]*X; + psF32 py = PAR[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; + + psF32 pr = PAR8*z; + psF32 p = pow(z, PAR[7] - 1.0); + psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); + psF32 f = PAR[1]*r + PAR[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 t = PAR[1]*r*r; + psF32 q = t*(PAR[7]*p + 4*PAR8*pr*pr*pr); + + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + deriv->data.F32[4] = -q*px*X; + deriv->data.F32[5] = -q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -t*log(z)*z*p; + } + return(f); +} + +psF64 psModelFlux_ZGAUSS(const psVector *params) +{ + float f, norm, z; + + psF32 *PAR = params->data.F32; + + psF64 A1 = PS_SQR(PAR[4]); + psF64 A2 = PS_SQR(PAR[5]); + psF64 A3 = PS_SQR(PAR[6]); + psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); + // Area is equivalent to 2 pi sigma^2 + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + psF32 pr = PAR8*z; + for (z = 0.005; z < 50; z += 0.01) { + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + norm += f; + } + norm *= 0.01; + + psF64 Flux = PAR[1] * Area * norm; + + return(Flux); +} + +// XXX need to define the radius along the major axis +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 psModelRadius_ZGAUSS (const psVector *params, psF64 flux) +{ + psF64 r, z, pr, f; + psF32 *PAR = params->data.F32; + + EllipseAxes axes; + EllipseShape shape; + + if (flux <= 0) return (1.0); + if (PAR[1] <= 0) return (1.0); + if (flux >= PAR[1]) return (1.0); + + // convert Sx,Sy,Sxy to major/minor axes + shape.sx = 1.0 / PAR[4]; + shape.sy = 1.0 / PAR[5]; + shape.sxy = PAR[6]; + + axes = EllipseShapeToAxes (shape); + psF64 dr = 1.0 / axes.major; + psF64 limit = flux / PAR[1]; + + // XXX : we can do this faster with an intelligent starting choice + for (r = 0.0; r < 20.0; r += dr) { + z = SQ(r); + pr = PAR8*z; + f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); + if (f < limit) break; + } + psF64 radius = 2.0 * axes.major * sqrt (z); + if (isnan(radius)) { + fprintf (stderr, "error in code\n"); + } + return (radius); +} + +bool psModelGuess_ZGAUSS (psModel *model, psSource *source) { + + psVector *params = model->params; + + EllipseAxes axes; + EllipseShape shape; + EllipseMoments moments; + + moments.x2 = PS_SQR(source->moments->Sx); + moments.y2 = PS_SQR(source->moments->Sy); + moments.xy = source->moments->Sxy; + + axes = EllipseMomentsToAxes(moments); + shape = EllipseAxesToShape(axes); + + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = 1.0 / shape.sx; + params->data.F32[5] = 1.0 / shape.sy; + params->data.F32[6] = shape.sxy; + params->data.F32[7] = 1.9; + return(true); +} + +bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) { + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + out[3] = in[3]; + + for (int i = 4; i < 8; i++) { + psPolynomial2D *poly = psf->params->data[i-4]; + out[i] = Polynomial2DEval (poly, out[2], out[3]); + } + return(true); +} Index: unk/psphot/src/onesource.c =================================================================== --- /trunk/psphot/src/onesource.c (revision 4953) +++ (revision ) @@ -1,43 +1,0 @@ -# include "psphot.h" - -bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y) { - - bool status; - int modelType; - float sky, obsMag, fitMag; - - modelType = psModelSetType (modelName); - if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName); - - float INNER = psMetadataLookupF32 (&status, config, "INNER_RADIUS"); - float OUTER = psMetadataLookupF32 (&status, config, "OUTER_RADIUS"); - float MRAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); - float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS"); - - psSource *source = pmSourceAlloc(); - source->peak = pmPeakAlloc (x, y, 0, 0); - pmSourceDefinePixels (source, imdata, x, y, OUTER); - - status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER); - status = pmSourceMoments_EAM (source, MRAD); - source->peak->counts = source->moments->Peak; - - psModel *model = pmSourceModelGuess (source, modelType); - fprintf (stderr, "guess slope: %f\n", model->params->data.F32[7]); - - psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80); - status = pmSourceFitModel_EAM (source, model, false); - - pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask); - fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag); - - // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - - DumpImage (source->pixels, "resid.fits"); - DumpImage (source->mask, "mask.fits"); - return true; -} Index: /trunk/psphot/src/pmModelInit.c =================================================================== --- /trunk/psphot/src/pmModelInit.c (revision 4954) +++ /trunk/psphot/src/pmModelInit.c (revision 4954) @@ -0,0 +1,111 @@ +# include +# include +# include "pmObjects_EAM.h" + +double hypot(double x, double y); +double sqrt (double x); + +# include "models/pmModel_GAUSS.c" +# include "models/pmModel_PGAUSS.c" +# include "models/pmModel_QGAUSS.c" +# include "models/pmModel_SGAUSS.c" + +static pmModelGroup models[] = { + {"PS_MODEL_GAUSS", 7, pmModelFunc_GAUSS, pmModelFlux_GAUSS, pmModelRadius_GAUSS, pmModelLimits_GAUSS, pmModelGuess_GAUSS, pmModelFromPSF_GAUSS, pmModelFitStatus_GAUSS}, + {"PS_MODEL_PGAUSS", 7, pmModelFunc_PGAUSS, pmModelFlux_PGAUSS, pmModelRadius_PGAUSS, pmModelLimits_PGAUSS, pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS}, + {"PS_MODEL_QGAUSS", 8, pmModelFunc_QGAUSS, pmModelFlux_QGAUSS, pmModelRadius_QGAUSS, pmModelLimits_QGAUSS, pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS}, + {"PS_MODEL_SGAUSS", 9, pmModelFunc_SGAUSS, pmModelFlux_SGAUSS, pmModelRadius_SGAUSS, pmModelLimits_SGAUSS, pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS}, +}; + +pmModelFunc pmModelFunc_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelFunc); +} + +pmModelFlux pmModelFlux_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelFlux); +} + +pmModelRadius pmModelRadius_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelRadius); +} + +pmModelLimits pmModelLimits_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelLimits); +} + +pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelGuessFunc); +} + +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelFitStatusFunc); +} + +pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].modelFromPSFFunc); +} + +psS32 pmModelParameterCount (pmModelType type) { + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (0); + } + return (models[type].nParams); +} + +psS32 pmModelSetType (char *name) { + + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + for (int i = 0; i < Nmodels; i++) { + if (!strcmp(models[i].name, name)) { + return (i); + } + } + return (-1); +} + +char *pmModelGetType (pmModelType type) { + + int Nmodels = sizeof (models) / sizeof (pmModelGroup); + if ((type < 0) || (type >= Nmodels)) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return (NULL); + } + return (models[type].name); +} Index: /trunk/psphot/src/pmObjects_EAM.c =================================================================== --- /trunk/psphot/src/pmObjects_EAM.c (revision 4953) +++ /trunk/psphot/src/pmObjects_EAM.c (revision 4954) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-06 08:05:08 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 03:33:01 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -66,6 +66,5 @@ pmModelAlloc(): Allocate the pmModel structure, along with its parameters, and initialize the type member. Initialize the params to 0.0. -XXX EAM: changing params and dparams to psVector -XXX EAM: simplifying code with psModelParameterCount +XXX EAM: simplifying code with pmModelParameterCount *****************************************************************************/ pmModel *pmModelAlloc(pmModelType type) @@ -76,7 +75,7 @@ tmp->chisq = 0.0; tmp->nIter = 0; - psS32 Nparams = psModelParameterCount (type); + psS32 Nparams = pmModelParameterCount (type); if (Nparams == 0) { - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); return(NULL); } @@ -167,5 +166,5 @@ // if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { + (vector->data.F32[0] > threshold)) { count++; } @@ -176,6 +175,6 @@ for (psU32 i = 1; i < n-1 ; i++) { if ((vector->data.F32[i] > vector->data.F32[i-1]) && - (vector->data.F32[i] > vector->data.F32[i+1]) && - (vector->data.F32[i] > threshold)) { + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { count++; } @@ -186,5 +185,5 @@ // if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { + (vector->data.F32[n-1] > threshold)) { count++; } @@ -201,5 +200,5 @@ // if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { + (vector->data.F32[0] > threshold)) { tmpVector->data.U32[count++] = 0; } @@ -210,6 +209,6 @@ for (psU32 i = 1; i < (n-1) ; i++) { if ((vector->data.F32[i] > vector->data.F32[i-1]) && - (vector->data.F32[i] > vector->data.F32[i+1]) && - (vector->data.F32[i] > threshold)) { + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { tmpVector->data.U32[count++] = i; } @@ -220,5 +219,5 @@ // if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { + (vector->data.F32[n-1] > threshold)) { tmpVector->data.U32[count++] = n-1; } @@ -319,6 +318,6 @@ if (col == 0) { if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { @@ -328,8 +327,8 @@ } else if (col < (image->numCols - 1)) { if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); @@ -339,6 +338,6 @@ } else if (col == (image->numCols - 1)) { if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { if (image->data.F32[row][col] > threshold) { list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); @@ -375,8 +374,8 @@ // If col==0, then we can not read col-1 pixels if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { myType = PM_PEAK_EDGE; list = myListAddPeak(list, row, col, image->data.F32[row][col], myType); @@ -385,31 +384,31 @@ // This is an interior pixel if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { myType = PM_PEAK_LONE; } if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || - (image->data.F32[row][col] == image->data.F32[row-1][col]) || - (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || - (image->data.F32[row][col] == image->data.F32[row][col-1]) || - (image->data.F32[row][col] == image->data.F32[row][col+1]) || - (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || - (image->data.F32[row][col] == image->data.F32[row+1][col]) || - (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { + (image->data.F32[row][col] == image->data.F32[row-1][col]) || + (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || + (image->data.F32[row][col] == image->data.F32[row][col-1]) || + (image->data.F32[row][col] == image->data.F32[row][col+1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col]) || + (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { myType = PM_PEAK_FLAT; } @@ -421,9 +420,9 @@ // If col==numCols - 1, then we can not read col+1 pixels if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col])) { + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col])) { myType = PM_PEAK_EDGE; list = myListAddPeak(list, row, col, image->data.F32[row][col], myType); @@ -448,6 +447,6 @@ if (col == 0) { if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1])) { + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1])) { if (image->data.F32[row][col] > threshold) { list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); @@ -456,8 +455,8 @@ } else if (col < (image->numCols - 1)) { if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1])) { + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1])) { if (image->data.F32[row][col] > threshold) { list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); @@ -467,6 +466,6 @@ } else if (col == (image->numCols - 1)) { if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row][col-1])) { + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1])) { if (image->data.F32[row][col] > threshold) { list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); @@ -495,5 +494,5 @@ (y >= valid->y0) && (y <= valid->y1)) { - return(true); + return(true); } @@ -511,4 +510,8 @@ XXX: changed API to create a NEW output psArray (should change name as well) + +XXX: Do we free the psList elements of those culled peaks? + +XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? *****************************************************************************/ psList *pmCullPeaks(psList *peaks, @@ -517,5 +520,4 @@ { PS_ASSERT_PTR_NON_NULL(peaks, NULL); - // PS_ASSERT_PTR_NON_NULL(valid, NULL); psListElem *tmpListElem = (psListElem *) peaks->head; @@ -526,6 +528,6 @@ pmPeak *tmpPeak = (pmPeak *) tmpListElem->data; if ((tmpPeak->counts > maxValue) || - ((valid != NULL) && - (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { + ((valid != NULL) && + (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { psListRemoveData(peaks, (psPtr) tmpPeak); } @@ -625,59 +627,19 @@ // XXX EAM : I added this code to stay on the image. So did George // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) + // XXX EAM : dropped the Annulus width (not needed) + // XXX EAM : dropped off-boundary tests (not needed) psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - outerRadiusS32); psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1); psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - outerRadiusS32); psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1); - // AnulusWidth == number of pixels width in the annulus. We add one since - // the pixels at the inner AND outher radius are included. - // XXX EAM : not used : psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32); - // Example: assume an outer/inner radius of 20/10. Then the subimage - // should have width/length of 40. An 18-by-18 interior region will - // be masked. - // printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32); - // printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32); - // printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth); - - // XXX EAM : these tests should not be needed: we can never hit this error because of above - # if (1) - - if (SubImageStartRow < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", - SubImageStartRow); - return(NULL); - } - if (SubImageEndRow > image->numRows) { - psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", - SubImageEndRow); - return(NULL); - } - if (SubImageStartCol < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", - SubImageStartCol); - return(NULL); - } - if (SubImageEndCol > image->numCols) { - psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", - SubImageEndCol); - return(NULL); - } - # endif - // // Grab a subimage of the original image of size (2 * outerRadius). // - // XXX: Must fix for new psImageSubset - // psImage *subImage = psImageSubset((psImage *) image, - // SubImageStartCol, - // SubImageStartRow, - // SubImageEndCol, - // SubImageEndRow); - // printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows); - psRegion tmpRegion = psRegionSet(SubImageStartCol, - SubImageEndCol, - SubImageStartRow, - SubImageEndRow); - psImage *subImage = psImageSubset((psImage *) image, tmpRegion); - // XXX EAM : can merge these: psImageSubset (image, psRegionSet ()); + // XXX EAM : merged psImageSubset & psRegionSet + // XXX EAM : cast for image is needed because of const (above) + psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, + SubImageEndCol, + SubImageStartRow, + SubImageEndRow)); psImage *subImageMask = psImageAlloc(subImage->numCols, @@ -685,14 +647,16 @@ PS_TYPE_U8); // XXX EAM : for consistency, mask needs col0,row0 set to match image - subImageMask->col0 = subImage->col0; - subImageMask->row0 = subImage->row0; - - // - // Loop through the subimage mask, initialize mask to 1 (invalid pixel) + // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!! + // subImageMask->col0 = subImage->col0; + // subImageMask->row0 = subImage->row0; + psAbort ("pmObjects", "must fix this error before psphot will work!!!\n"); + + // + // Loop through the subimage mask, initialize mask to 0x01 (invalid pixel) // XXX EAM : use PSPHOT_MASK_INVALID? // for (psS32 row = 0 ; row < subImageMask->numRows; row++) { for (psS32 col = 0 ; col < subImageMask->numCols; col++) { - subImageMask->data.U8[row][col] = 1; + subImageMask->data.U8[row][col] = 0x01; } } @@ -702,4 +666,6 @@ // XXX this uses a static mask value of 0 // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image + // XXX EAM : this masks the pixels in the center region + // XXX EAM : should we be using psImageMaskRegion ??? psS32 StartRow = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32); psS32 EndRow = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1); @@ -712,5 +678,5 @@ } - + // XXX EAM : make this trace information?? // for (psS32 row = 0 ; row < subImage->numRows; row++) { // for (psS32 col = 0 ; col < subImage->numCols; col++) { @@ -800,5 +766,5 @@ XXX: mask values? *****************************************************************************/ -bool pmSourceMoments(psSource *source, +bool pmSourceMoments(pmSource *source, psF32 radius) { @@ -849,5 +815,5 @@ psS32 imgColCoord = col + source->pixels->col0; psS32 imgRowCoord = row + source->pixels->row0; - if (CheckRadius(source->peak, + if (checkRadius(source->peak, radius, imgColCoord, @@ -874,6 +840,6 @@ // XXX EAM - the limit is a bit arbitrary. make it user defined? if ((numPixels < 3) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); - return (false); + psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); + return (false); } @@ -890,8 +856,8 @@ y = Y1/Sum; if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace (".psModules.pmSourceMoments", 5, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - return (false); + psTrace (".psModules.pmSourceMoments", 5, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + return (false); } @@ -914,39 +880,4 @@ return(true); - - // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it - # if (0) - // - // second loop: get the difference sums - // - X2 = Y2 = 0; - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { - psS32 imgColCoord = col + source->pixels->col0; - psS32 imgRowCoord = row + source->pixels->row0; - if (checkRadius(source->peak, - radius, - imgColCoord, - imgRowCoord)) { - psF32 xDiff = (psF32) (imgColCoord - source->peak->x); - psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - - Sum+= pDiff; - X2+= PS_SQR(xDiff - x) * pDiff; - Y2+= PS_SQR(yDiff - y) * pDiff; - } - } - } - } - - // - // second moment X = sqrt (X2/Sum) - // - source->moments->Sx = (X2/Sum); - source->moments->Sy = (Y2/Sum); - return(true); - # endif } @@ -990,6 +921,6 @@ { - # define NPIX 10 - # define SCALE 0.1 +# define NPIX 10 +# define SCALE 0.1 psArray *peaks = NULL; @@ -997,6 +928,6 @@ pmPSFClump psfClump = emptyClump; - // PS_PTR_CHECK_NULL(sources, emptyClump); - // PS_PTR_CHECK_NULL(metadata, emptyClump); + PS_ASSERT_PTR_NON_NULL(sources, emptyClump); + PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); // find the sigmaX, sigmaY clump @@ -1010,5 +941,5 @@ splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); for (int i = 0; i < splane->numRows; i++) { - memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); + memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); } @@ -1139,4 +1070,6 @@ XXX: How can this function ever return FALSE? + +XXX EAM : add the saturated mask value to metadata *****************************************************************************/ @@ -1178,10 +1111,11 @@ psF32 A = 4 * M_PI * sigX * sigY; psF32 B = tmpSrc->moments->Sky; - psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN))); - psF32 SN = (S * PS_SQRT_F32(GAIN) / RT); + psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN))); + psF32 SN = (S * sqrt(GAIN) / RT); tmpSrc->moments->SN = SN; - // XXX EAM : can we use the value of SATIRATE if mask is NULL? - int Nsatpix = pmCountSatPixels (tmpSrc->mask); + // XXX EAM : can we use the value of SATURATE if mask is NULL? +# define MASK_SATURATED 0x02 + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, MASK_SATURATED); // saturated star (size consistent with PSF or larger) @@ -1253,18 +1187,4 @@ } -int pmCountSatPixels (psImage *image) -{ - int Nsatpix = 0; - - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - if (image->data.U8[i][j] & 0x02) { - Nsatpix ++; - } - } - } - return (Nsatpix); -} - /****************************************************************************** pmSourceSetPixelsCircle(source, image, radius) @@ -1286,5 +1206,5 @@ PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(source->moments, false); - // PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_FLOAT_COMPARE(0.0, radius, false); @@ -1305,29 +1225,4 @@ psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); - // XXX EAM : this should not be needed: we can never hit this error - # if (1) - - if (SubImageStartRow < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", - SubImageStartRow); - return(false); - } - if (SubImageEndRow > image->numRows) { - psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", - SubImageEndRow); - return(false); - } - if (SubImageStartCol < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", - SubImageStartCol); - return(false); - } - if (SubImageEndCol > image->numCols) { - psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", - SubImageEndCol); - return(false); - } - # endif - // XXX: Must recycle image. // XXX EAM: this message reflects a programming error we know about. @@ -1338,10 +1233,4 @@ psFree(source->pixels); } - // XXX: Must fix this. psImageSubset() has different parameters in latest CVS. - // source->pixels = psImageSubset((psImage *) image, - // SubImageStartCol, - // SubImageStartRow, - // SubImageEndCol, - // SubImageEndRow); source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, SubImageStartRow, @@ -1390,14 +1279,13 @@ functions will be in image, not subImage coords. Remember this. *****************************************************************************/ -psModel *pmSourceModelGuess(psSource *source, - psModelType modelType) -{ - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - - psModel *model = psModelAlloc(modelType); - - psModelGuessFunc modelGuessFunc = psModelGuessFunc_GetFunction (modelType); +pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) +{ + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction (modelType); modelGuessFunc (model, source); return(model); @@ -1425,21 +1313,13 @@ testing. Try to reproduce that and debug. *****************************************************************************/ -static psF32 evalModel(pmSource *src, - psU32 row, - psU32 col) -{ - PS_ASSERT_PTR_NON_NULL(src, false); - PS_ASSERT_PTR_NON_NULL(src->modelPSF, false); - PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false); - -XXX EAM : I've made this a public function -XXX EAM : this now uses a psModel as the input -XXX EAM : it was using src->type to find the model, not model->type -*****************************************************************************/ -psF32 psModelEval(psModel *model, psImage *image, psS32 col, psS32 row) -{ - PS_PTR_CHECK_NULL(image, false); - PS_PTR_CHECK_NULL(model, false); - PS_PTR_CHECK_NULL(model->params, false); + +// XXX EAM : I have made this a public function +// XXX EAM : this now uses a pmModel as the input +// XXX EAM : it was using src->type to find the model, not model->type +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) +{ + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_PTR_NON_NULL(model->params, false); // Allocate the x coordinate structure and convert row/col to image space. @@ -1449,7 +1329,7 @@ x->data.F32[1] = (psF32) (row + image->row0); psF32 tmpF; - psModelFunc modelFunc; - - modelFunc = psModelFunc_GetFunction (model->type); + pmModelFunc modelFunc; + + modelFunc = pmModelFunc_GetFunction (model->type); tmpF = modelFunc (NULL, model->params, x); psFree(x); @@ -1491,7 +1371,7 @@ } - // XXX EAM : i changed this to match psModelEval above, but see + // XXX EAM : i changed this to match pmModelEval above, but see // XXX EAM the note below in pmSourceContour - psF32 oldValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow); + psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { return(((psF32) (subCol + source->pixels->col0))); @@ -1514,5 +1394,5 @@ while (subCol != lastColumn) { - psF32 newValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow); + psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { return((psF32) (subCol + source->pixels->col0)); @@ -1674,5 +1554,5 @@ // tests below could be conditions (!NULL) - psModelFunc modelFunc = psModelFunc_GetFunction (model->type); + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); psVector *params = model->params; @@ -1694,6 +1574,6 @@ } if (count < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - return(false); + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + return(false); } @@ -1722,15 +1602,15 @@ psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); + PM_SOURCE_FIT_MODEL_TOLERANCE); // PSF model only fits first 4 parameters, FLT model fits all if (PSF) { - paramMask = psVectorAlloc (params->n, PS_TYPE_U8); - for (int i = 0; i < 4; i++) { - paramMask->data.U8[i] = 0; - } - for (int i = 4; i < paramMask->n; i++) { - paramMask->data.U8[i] = 1; - } + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } } @@ -1782,17 +1662,16 @@ bool p_pmSourceAddOrSubModel(psImage *image, psImage *mask, - psModel *model, + pmModel *model, bool center, psS32 flag) { - // XXX EAM : convert to ASSERTS - // PS_PTR_CHECK_NULL(model, false); - // PS_IMAGE_CHECK_NULL(image, false); - // PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); psVector *x = psVectorAlloc(2, PS_TYPE_F32); psVector *params = model->params; - psModelFunc modelFunc = psModelFunc_GetFunction (model->type); + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); psS32 imageCol; psS32 imageRow; @@ -1842,5 +1721,5 @@ bool pmSourceAddModel(psImage *image, psImage *mask, - psModel *model, + pmModel *model, bool center) { @@ -1852,5 +1731,5 @@ bool pmSourceSubModel(psImage *image, psImage *mask, - psModel *model, + pmModel *model, bool center) { Index: /trunk/psphot/src/pmObjects_EAM.h =================================================================== --- /trunk/psphot/src/pmObjects_EAM.h (revision 4953) +++ /trunk/psphot/src/pmObjects_EAM.h (revision 4954) @@ -1,3 +1,101 @@ - +/** @file pmObjects.h + * + * This file will ... + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-07 03:33:01 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_OBJECTS_H) +#define PM_OBJECTS_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" + +// XXX EAM : this function should be added to psLib +int psImageCountPixelMask (psImage *mask, psU8 value); + +// XXX EAM : bug : missing from psLib/*.h +bool psMinimizeGaussNewtonDelta (psVector *delta, + const psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func); + +// XXX EAM : switch to psPoly when ready +psF64 Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y); + +// XXX EAM : psEllipse needs to be be included in psLib +# include "psEllipse.h" + +/** pmPeakType + * + * A peak pixel may have several features which may be determined when the + * peak is found or measured. These are specified by the pmPeakType enum. + * PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate + * neighbors. The PM_PEAK_EDGE represents a peak pixel which touching the image + * edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific + * number of neighbors at the same value, within some tolarence: + * + */ +typedef enum { + PM_PEAK_LONE, ///< Isolated peak. + PM_PEAK_EDGE, ///< Peak on edge. + PM_PEAK_FLAT, ///< Peak has equal-value neighbors. + PM_PEAK_UNDEF ///< Undefined. +} pmPeakType; + +/** pmPeak data structure + * + * + * + */ +typedef struct +{ + int x; ///< X-coordinate of peak pixel. + int y; ///< Y-coordinate of peak pixel. + float counts; ///< Value of peak pixel (above sky?). + pmPeakType class; ///< Description of peak. +} +pmPeak; + +/** pmMoments data structure + * + * + * + */ +typedef struct +{ + float x; ///< X-coord of centroid. + float y; ///< Y-coord of centroid. + float Sx; ///< x-second moment. + float Sy; ///< y-second moment. + float Sxy; ///< xy cross moment. + float Sum; ///< Pixel sum above sky (background). + float Peak; ///< Peak counts above sky. + float Sky; ///< Sky level (background). + float SN; ///< approx signal-to-noise + int nPixels; ///< Number of pixels used. +} pmMoments; + +typedef int pmModelType; + +/** pmPSFClump data structure + * + * + * + */ typedef struct { @@ -8,33 +106,26 @@ } pmPSFClump; -typedef struct -{ - float x; ///< X-coord of centroid. - float y; ///< Y-coord of centroid. - float Sx; ///< x-second moment. - float Sy; ///< y-second moment. - float Sxy; ///< xy cross moment. - float Sum; ///< Pixel sum above sky (background). - float Peak; ///< Peak counts above sky. - float Sky; ///< Sky level (background). - float SN; - int nPixels; ///< Number of pixels used. -} -pmMoments; - -typedef psS32 pmModelType; - -typedef struct -{ - pmModelType type; ///< Model to be used. - psVector *params; ///< Paramater values. - psVector *dparams; ///< Parameter errors. - float radius; ///< fit radius actually used - float chisq; ///< Fit chi-squared. - int nDOF; ///< number of degrees of freedom - int nIter; ///< number of iterations to reach min +/** pmModel data structure + * + * + * + */ +typedef struct +{ + pmModelType type; ///< Model to be used. + psVector *params; ///< Paramater values. + psVector *dparams; ///< Parameter errors. + float chisq; ///< Fit chi-squared. + int nDOF; ///< number of degrees of freedom + int nIter; ///< number of iterations to reach min + float radius; ///< fit radius actually used } pmModel; +/** pmSourceType enumeration + * + * + * + */ typedef enum { PS_SOURCE_DEFECT, @@ -58,97 +149,267 @@ } pmSourceType; -typedef struct -{ - pmPeak *peak; ///< Description of peak pixel. - psImage *pixels; ///< Rectangular region including object pixels. - psImage *noise; ///< Mask which marks pixels associated with objects. - psImage *mask; ///< Mask which marks pixels associated with objects. - pmMoments *moments; ///< Basic moments measure for the object. - pmModel *modelPSF; ///< PSF Model fit (parameters and type) - pmModel *modelFLT; ///< FLT Model fit (parameters and type). - pmSourceType type; ///< Best identification of object. +/** pmSource data structure + * + * This source has the capacity for several types of measurements. The + * simplest measurement of a source is the location and flux of the peak pixel + * associated with the source: + * + */ +typedef struct +{ + pmPeak *peak; ///< Description of peak pixel. + psImage *pixels; ///< Rectangular region including object pixels. + psImage *noise; ///< Mask which marks pixels associated with objects. + psImage *mask; ///< Mask which marks pixels associated with objects. + pmMoments *moments; ///< Basic moments measure for the object. + pmModel *modelPSF; ///< PSF Model fit (parameters and type) + pmModel *modelFLT; ///< FLT Model fit (parameters and type). + pmSourceType type; ///< Best identification of object. } pmSource; -// XXX EAM function to return pmModel -typedef psMinimizeLMChi2Func pmModelFunc; -typedef psF64 (*pmModelFlux)(const psVector *params); -typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); -typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); -typedef psF64 (*pmModelRadius)(const psVector *params, double flux); - -typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - -// XXX EAM structure to carry model options -typedef struct { - char *name; - int nParams; - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuessFunc; - pmModelFromPSFFunc modelFromPSFFunc; - pmModelFitStatusFunc modelFitStatusFunc; -} pmModelGroup; - -// XXX EAM : strucures to define elliptical shape parameters -typedef struct { - double major; - double minor; - double theta; -} EllipseAxes; - -typedef struct { - double x2; - double y2; - double xy; -} EllipseMoments; - -typedef struct { - double sx; - double sy; - double sxy; -} EllipseShape; - -EllipseAxes EllipseMomentsToAxes (EllipseMoments moments); -EllipseShape EllipseAxesToShape (EllipseAxes axes); -EllipseAxes EllipseShapeToAxes (EllipseShape shape); +/** pmPSF data structure + * + * + * + */ +typedef struct +{ + pmModelType type; ///< PSF Model in use + psArray *params; ///< Model parameters (psPolynomial2D) + float chisq; ///< PSF goodness statistic + int nPSFstars; ///< number of stars used to measure PSF +} +pmPSF; + +pmPeak *pmPeakAlloc( + int x, ///< Row-coordinate in image space + int y, ///< Col-coordinate in image space + float counts, ///< The value of the peak pixel + pmPeakType class ///< The type of peak pixel +); + +pmMoments *pmMomentsAlloc(); +pmModel *pmModelAlloc(pmModelType type); +pmSource *pmSourceAlloc(); + +/****************************************************************************** +pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector +above the given threshold. Returns a vector of type PS_TYPE_U32 containing +the location (x value) of all peaks. +*****************************************************************************/ +psVector *pmFindVectorPeaks( + const psVector *vector, ///< The input vector (float) + float threshold ///< Threshold above which to find a peak +); + +/****************************************************************************** +pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage +above the given threshold. Returns a psList containing the location (x/y +value) of all peaks. +*****************************************************************************/ +psArray *pmFindImagePeaks( + const psImage *image, ///< The input image where peaks will be found (float) + float threshold ///< Threshold above which to find a peak +); + +/****************************************************************************** +pmCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have +a peak value above the given maximum, or fall outside the valid region. +*****************************************************************************/ +psList *pmCullPeaks( + psList *peaks, ///< The psList of peaks to be culled + float maxValue, ///< Cull peaks above this value + const psRegion *valid ///< Cull peaks otside this psRegion +); + +/****************************************************************************** +pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): + +*****************************************************************************/ +pmSource *pmSourceLocalSky( + const psImage *image, ///< The input image (float) + const pmPeak *peak, ///< The peak for which the pmSource struct is created. + psStatsOptions statsOptions, ///< The statistic used in calculating the background sky + float innerRadius, ///< The inner radius of the suqare annulus for calculating sky + float outerRadius ///< The outer radius of the suqare annulus for calculating sky +); + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceMoments( + pmSource *source, ///< The input pmSource for which moments will be computed + float radius ///< Use a circle of pixels around the peak +); /****************************************************************************** pmSourcePSFClump(pmArray *source, psMetaDeta *metadata): find the source PSF clump - *****************************************************************************/ -pmPSFClump pmSourcePSFClump(psArray *source, ///< The input psSource - psMetadata *metadata ///< Contains classification parameters - ); +*****************************************************************************/ +pmPSFClump pmSourcePSFClump( + psArray *source, ///< The input pmSource + psMetadata *metadata ///< Contains classification parameters +); /****************************************************************************** pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the source classification. - *****************************************************************************/ -bool pmSourceRoughClass(psArray *source, ///< The input psSource - psMetadata *metadata, ///< Contains classification parameters - pmPSFClump clump ///< Statistics about the PSF clump - ); - -pmModel *pmSourceModelGuess(pmSource *source, ///< The input psSource - pmModelType model ///< The type of model to be created. - ); - -// XXX EAM : added utility functions -pmModelType pmModelSetType (char *name); -pmModelFunc pmModelFunc_GetFunction (pmModelType type); -pmModelFlux pmModelFlux_GetFunction (pmModelType type); -pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type); -pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type); -pmModelRadius pmModelRadius_GetFunction (pmModelType type); -char *pmModelGetType (pmModelType type); -psS32 pmModelParameterCount (pmModelType type); -pmModelLimits pmModelLimits_GetFunction (pmModelType type); -pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); -int pmCountSatPixels (psImage *image); - -pmMoments *pmMomentsAlloc(); -pmModel *pmModelAlloc(pmModelType type); -pmSource *pmSourceAlloc(); +*****************************************************************************/ +bool pmSourceRoughClass( + psArray *source, ///< The input pmSource + psMetadata *metadata, ///< Contains classification parameters + pmPSFClump clump ///< Statistics about the PSF clump +); + +/****************************************************************************** +pmSourceSetPixelCircle(source, image, radius) +*****************************************************************************/ +bool pmSourceSetPixelsCircle( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) + float radius ///< The radius of the circle +); + +/****************************************************************************** + *****************************************************************************/ +pmModel *pmSourceModelGuess( + pmSource *source, ///< The input pmSource + pmModelType model ///< The type of model to be created. +); + +/****************************************************************************** + *****************************************************************************/ +typedef enum { + PS_CONTOUR_CRUDE, +} pmContourType; + +psArray *pmSourceContour( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) (this arg should be removed) + float level, ///< The level of the contour + pmContourType mode ///< Currently this must be PS_CONTOUR_CRUDE +); + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceFitModel( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceAddModel( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center ///< A boolean flag that determines whether pixels are centered +); + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceSubModel( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center ///< A boolean flag that determines whether pixels are centered +); + +/** + * + * The object model functions are defined to allow for the flexible addition + * of new object models. Every object model, with parameters represented by + * pmModel, has an associated set of functions which provide necessary support + * operations. A set of abstract functions allow the programmer to select the + * approriate function or property for a specific named object model. + * + */ + +/** + * + * This function is the model chi-square minimization function for this model. + * + */ +typedef psMinimizeLMChi2Func pmModelFunc; + + +/** + * + * This function returns the integrated flux for the given model parameters. + */ +typedef psF64 (*pmModelFlux)(const psVector *params); + + +/** + * + * This function returns the radius at which the given model and parameters + * achieves the given flux. + * + */ +typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + +/** + * + * This function sets the model parameter limits vectors for the given model + * + */ +typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); + +/** + * + * This function provides the model guess parameters based on the details of + * the given source. + * + */ +typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + + +/** + * + * This function constructs the PSF model for the given source based on the + * supplied psf and the FLT model for the object. + * + */ +typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); + + +/** + * + * This function returns the success / failure status of the given model fit + * + */ +typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +/** + * + * Each of the function types above has a corresponding function which returns + * the function given the model type: + * + */ +pmModelFunc pmModelFunc_GetFunction (pmModelType type); +pmModelFlux pmModelFlux_GetFunction (pmModelType type); +pmModelRadius pmModelRadius_GetFunction (pmModelType type); +pmModelLimits pmModelLimits_GetFunction (pmModelType type); +pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type); +pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type); +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); + +// pmModelGroup utility functions +int pmModelParameterCount (pmModelType type); +char *pmModelGetType (pmModelType type); +pmModelType pmModelSetType (char *name); + +// structure to carry model group functions +typedef struct { + char *name; + int nParams; + pmModelFunc modelFunc; + pmModelFlux modelFlux; + pmModelRadius modelRadius; + pmModelLimits modelLimits; + pmModelGuessFunc modelGuessFunc; + pmModelFromPSFFunc modelFromPSFFunc; + pmModelFitStatusFunc modelFitStatusFunc; +} pmModelGroup; + +#endif Index: /trunk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 4953) +++ /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 4954) @@ -16,5 +16,5 @@ psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); - psImageSmooth (smooth, SIGMA, NSIGMA); + psImageSmooth_EAM (smooth, SIGMA, NSIGMA); psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); Index: /trunk/psphot/src/pmSourceUtils.c =================================================================== --- /trunk/psphot/src/pmSourceUtils.c (revision 4953) +++ /trunk/psphot/src/pmSourceUtils.c (revision 4954) @@ -7,20 +7,14 @@ psF32 Radius) { - psRegion *srcRegion; + psRegion srcRegion; // Grab a subimage of the original image of size (2 * outerRadius). - srcRegion = psRegionSquare (x, y, Radius); - srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion); + srcRegion = psRegionForSquare (x, y, Radius); + srcRegion = psRegionForImage (imdata->image, &srcRegion); - // use these when psImageSubset is updated - // psImage *subImage = psImageSubset(imdata->image, srcRegion); - // psImage *subImageMask = psImageSubset(imdata->mask, srcRegion); - // psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion); + mySource->pixels = psImageSubset(imdata->image, srcRegion); + mySource->noise = psImageSubset(imdata->noise, srcRegion); + mySource->mask = psImageSubset(imdata->mask, srcRegion); - mySource->pixels = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - mySource->noise = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - mySource->mask = psImageSubset(imdata->mask, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - - psFree (srcRegion); return(mySource); } Index: /trunk/psphot/src/psEllipse.c =================================================================== --- /trunk/psphot/src/psEllipse.c (revision 4954) +++ /trunk/psphot/src/psEllipse.c (revision 4954) @@ -0,0 +1,59 @@ +# include "pslib.h" +# include "psEllipse.h" + +EllipseAxes EllipseMomentsToAxes (EllipseMoments moments) { + + EllipseAxes axes; + + double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy)); + if (f > (moments.x2 + moments.y2) / 2.0) { + f = 0.98*(moments.x2 + moments.y2) / 2.0; + } + + axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f); + axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f); + axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2; + // theta in radians + + return (axes); +} + +EllipseShape EllipseAxesToShape (EllipseAxes axes) { + + EllipseShape shape; + + double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor); + double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor); + + double sxr = r1 + r2*cos(2*axes.theta); + double syr = r1 - r2*cos(2*axes.theta); + + shape.sx = 1.0 / sqrt(sxr); + shape.sy = 1.0 / sqrt(syr); + shape.sxy = r2*sin(2*axes.theta); + + return (shape); +} + +EllipseAxes EllipseShapeToAxes (EllipseShape shape) { + + EllipseAxes axes; + + double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy); + double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy); + + // force the axis ratio to be less than 10 + double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2)); + + shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1); + + axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0; + + double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); + double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); + + axes.minor = 1.0 / sqrt (Ar); + axes.major = 1.0 / sqrt (Br); + + return (axes); +} Index: /trunk/psphot/src/psEllipse.h =================================================================== --- /trunk/psphot/src/psEllipse.h (revision 4954) +++ /trunk/psphot/src/psEllipse.h (revision 4954) @@ -0,0 +1,24 @@ +// strucures to define elliptical shape parameters +typedef struct { + double major; + double minor; + double theta; +} EllipseAxes; + +typedef struct { + double x2; + double y2; + double xy; +} EllipseMoments; + +typedef struct { + double sx; + double sy; + double sxy; +} EllipseShape; + +// conversions between elliptical shape representations +EllipseAxes EllipseMomentsToAxes (EllipseMoments moments); +EllipseShape EllipseAxesToShape (EllipseAxes axes); +EllipseAxes EllipseShapeToAxes (EllipseShape shape); + Index: /trunk/psphot/src/psImageData.c =================================================================== --- /trunk/psphot/src/psImageData.c (revision 4953) +++ /trunk/psphot/src/psImageData.c (revision 4954) @@ -20,5 +20,5 @@ imdata->mask = mask; - p_psMemSetDeallocator(imdata, (psFreeFcn) psImageDataFree); + psMemSetDeallocator(imdata, (psFreeFunc) psImageDataFree); return (imdata); } Index: /trunk/psphot/src/psLibUtils.c =================================================================== --- /trunk/psphot/src/psLibUtils.c (revision 4953) +++ /trunk/psphot/src/psLibUtils.c (revision 4954) @@ -1,3 +1,7 @@ -# include "psphot.h" +# include // for strncasecmp +# include + +// 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 @@ -29,5 +33,5 @@ if (timers == NULL) timers = psHashAlloc (16); - start = psTimeGetTime (PS_TIME_UTC); + start = psTimeGetNow (PS_TIME_UTC); psHashAdd (timers, name, start); psFree (start); @@ -47,5 +51,5 @@ if (start == NULL) return (0); - mark = psTimeGetTime (PS_TIME_UTC); + mark = psTimeGetNow (PS_TIME_UTC); delta = psTimeDelta (mark, start); psFree (mark); @@ -53,4 +57,32 @@ return (delta); +} + +// 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); +} + +// 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); } @@ -114,38 +146,11 @@ } -// 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); -} - -// 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); -} - // alternate implementation of this function from pmObjects.c +// XXX EAM : move this to pmObjects_EAM.c psVector *psGetRowVectorFromImage(psImage *image, psU32 row) { - PS_IMAGE_CHECK_NULL(image, NULL); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); @@ -201,5 +206,6 @@ } -void psImageSmooth (psImage *image, float sigma, float Nsigma) { +// XXX EAM : this is now in psLib +void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma) { int Nx, Ny, Npixel, Nrange; @@ -314,5 +320,22 @@ } +// 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; @@ -321,4 +344,5 @@ return (region); } +# endif // set actual region based on image parameters: @@ -327,7 +351,9 @@ // 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, us the child Nx,Ny -// force range to be on this subimage -// XXX EAM : this needs to be changes to use psRegion rather than psRegion* +// 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) { @@ -356,7 +382,10 @@ 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) { @@ -375,7 +404,10 @@ } } +# 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) { @@ -396,7 +428,10 @@ } } +# 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) { @@ -419,7 +454,10 @@ } } +# 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) { @@ -442,4 +480,5 @@ } } +# endif psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { @@ -455,2 +494,19 @@ 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: /trunk/psphot/src/psLine.c =================================================================== --- /trunk/psphot/src/psLine.c (revision 4953) +++ /trunk/psphot/src/psLine.c (revision 4954) @@ -14,5 +14,5 @@ psLine *line; line = psAlloc (sizeof(psLine)); - p_psMemSetDeallocator(line, (psFreeFcn) psLineFree); + psMemSetDeallocator(line, (psFreeFunc) psLineFree); line->Nline = 0; Index: /trunk/psphot/src/psMinimize.c =================================================================== --- /trunk/psphot/src/psMinimize.c (revision 4953) +++ /trunk/psphot/src/psMinimize.c (revision 4954) @@ -13,12 +13,12 @@ psMinimizeLMChi2Func func) { - PS_PTR_CHECK_NULL(min, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_PTR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_NULL(y, NULL); - PS_VECTOR_CHECK_EMPTY(y, NULL); - PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); - PS_PTR_CHECK_NULL(func, NULL); + // PS_PTR_CHECK_NULL(min, NULL); + // PS_VECTOR_CHECK_NULL(params, NULL); + // PS_VECTOR_CHECK_EMPTY(params, NULL); + // PS_PTR_CHECK_NULL(x, NULL); + // PS_VECTOR_CHECK_NULL(y, NULL); + // PS_VECTOR_CHECK_EMPTY(y, NULL); + // PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); + // PS_PTR_CHECK_NULL(func, NULL); // this function has test and current values for several things Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 4953) +++ /trunk/psphot/src/psModulesUtils.c (revision 4954) @@ -11,20 +11,20 @@ psImage *image = source->pixels; psImage *mask = source->mask; - psPeak *peak = source->peak; - psRegion *srcRegion; - - srcRegion = psRegionSquare (peak->x, peak->y, Radius); - srcRegion = psRegionForImage (srcRegion, mask, srcRegion); - psImageMaskRegion (mask, srcRegion, OR, 0x80); - + pmPeak *peak = source->peak; + psRegion srcRegion; + + // XXX EAM : psRegionXXX funcs need value not ptr + + srcRegion = psRegionForSquare (peak->x, peak->y, Radius); + srcRegion = psRegionForImage (mask, &srcRegion); + + psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED); psStats *myStats = psStatsAlloc(statsOptions); myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion (mask, srcRegion, AND, 0x7f); + psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED); psF64 tmpF64; p_psGetStatValue(myStats, &tmpF64); - psFree (myStats); - psFree (srcRegion); if (isnan(tmpF64)) return (false); @@ -40,10 +40,10 @@ const bool PSF) { - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_PTR_CHECK_NULL(source->pixels, false); - PS_PTR_CHECK_NULL(source->mask, false); - PS_PTR_CHECK_NULL(source->noise, false); + // PS_PTR_CHECK_NULL(source, false); + // PS_PTR_CHECK_NULL(source->moments, false); + // PS_PTR_CHECK_NULL(source->peak, false); + // PS_PTR_CHECK_NULL(source->pixels, false); + // PS_PTR_CHECK_NULL(source->mask, false); + // PS_PTR_CHECK_NULL(source->noise, false); psBool fitStatus = true; psBool onPic = true; @@ -188,4 +188,6 @@ XXX: mask values? + +XXX EAM : this version clips input pixels on S/N *****************************************************************************/ # define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) @@ -194,8 +196,8 @@ psF32 radius) { - PS_PTR_CHECK_NULL(source, NULL); - PS_PTR_CHECK_NULL(source->peak, NULL); - PS_PTR_CHECK_NULL(source->pixels, NULL); - PS_PTR_CHECK_NULL(source->noise, NULL); + // PS_PTR_CHECK_NULL(source, NULL); + // PS_PTR_CHECK_NULL(source->peak, NULL); + // PS_PTR_CHECK_NULL(source->pixels, NULL); + // PS_PTR_CHECK_NULL(source->noise, NULL); PS_FLOAT_COMPARE(0.0, radius, NULL); Index: /trunk/psphot/src/psPolynomials.c =================================================================== --- /trunk/psphot/src/psPolynomials.c (revision 4953) +++ /trunk/psphot/src/psPolynomials.c (revision 4954) @@ -1,3 +1,10 @@ # include "psphot.h" + +// XXX EAM : this file defines alternate versions of the polynomial fitting +// functions. these differ from psLib v.0.5.0 in the definition +// of nOrder (here we use nOrder in the mathematical sense, not +// like IRAF. ie, a linear function is of nOrder 1) +// psLib v 0.7.0 should be in sync with out defintion of nOrder +// If so, we can drop these functions and use the psLib versions. // write out the terms of the given 1D polynomial @@ -9,9 +16,9 @@ } -psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly) +psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly) { psS32 loop_x = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) { @@ -28,7 +35,7 @@ const psVector *x) { - PS_POLY_CHECK_NULL(myPoly, NULL); - PS_VECTOR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); + // PS_POLY_CHECK_NULL(myPoly, NULL); + // PS_VECTOR_CHECK_NULL(x, NULL); + // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); psVector *tmp; @@ -56,5 +63,5 @@ psPolynomialType type) { - PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL); + // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL); psS32 i = 0; @@ -65,6 +72,6 @@ newPoly->type = type; newPoly->n = nOrder; - newPoly->coeff = (psF32 *)psAlloc(nTerm * sizeof(psF32)); - newPoly->coeffErr = (psF32 *)psAlloc(nTerm * sizeof(psF32)); + newPoly->coeff = (psF64 *)psAlloc(nTerm * sizeof(psF64)); + newPoly->coeffErr = (psF64 *)psAlloc(nTerm * sizeof(psF64)); newPoly->mask = (psU8 *)psAlloc(nTerm * sizeof(psU8)); for (i = 0; i < nTerm; i++) { @@ -73,5 +80,5 @@ newPoly->mask[i] = 0; } - p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial1DFree); + psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial1DFree); return(newPoly); } @@ -120,4 +127,6 @@ // dump minutiae # ifndef PS_NO_TRACE +# if (0) + // XXX EAM : change from FILE to fd breaks this code: if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) { FILE *f = psTraceGetDestination (); @@ -127,4 +136,5 @@ } } +# endif # endif @@ -214,9 +224,9 @@ { - PS_POLY_CHECK_NULL(myPoly, NULL); - PS_VECTOR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(y, NULL); - PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL); + // PS_POLY_CHECK_NULL(myPoly, NULL); + // PS_VECTOR_CHECK_NULL(x, NULL); + // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL); + // PS_VECTOR_CHECK_NULL(y, NULL); + // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL); psVector *tmp; @@ -246,9 +256,9 @@ { - PS_POLY_CHECK_NULL(myPoly, NULL); - PS_VECTOR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); - PS_VECTOR_CHECK_NULL(y, NULL); - PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL); + // PS_POLY_CHECK_NULL(myPoly, NULL); + // PS_VECTOR_CHECK_NULL(x, NULL); + // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL); + // PS_VECTOR_CHECK_NULL(y, NULL); + // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL); psVector *tmp; @@ -291,6 +301,6 @@ psPolynomialType type) { - PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL); - PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL); + // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL); + // PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL); psS32 x = 0; @@ -305,11 +315,13 @@ newPoly->nY = nYorder; - newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *)); - newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *)); - newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *)); + newPoly->coeff = (psF64 **)psAlloc(nXterm * sizeof(psF64 *)); + newPoly->coeffErr = (psF64 **)psAlloc(nXterm * sizeof(psF64 *)); + newPoly->mask = (char **)psAlloc(nXterm * sizeof(char *)); + // XXX EAM : this is an error in the definition of the polynomial mask + // (should be psU8 not char) for (x = 0; x < nXterm; x++) { - newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32)); - newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32)); - newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8)); + newPoly->coeff[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64)); + newPoly->coeffErr[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64)); + newPoly->mask[x] = (char *)psAlloc(nYterm * sizeof(char)); } for (x = 0; x < nXterm; x++) { @@ -320,5 +332,5 @@ } } - p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial2DFree); + psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial2DFree); return(newPoly); } @@ -510,9 +522,9 @@ const psVector* dz) { - PS_VECTOR_CHECK_NULL(mask, NULL); - PS_VECTOR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_NULL(y, NULL); - PS_VECTOR_CHECK_NULL(z, NULL); - PS_VECTOR_CHECK_NULL(dz, NULL); + // PS_VECTOR_CHECK_NULL(mask, NULL); + // PS_VECTOR_CHECK_NULL(x, NULL); + // PS_VECTOR_CHECK_NULL(y, NULL); + // PS_VECTOR_CHECK_NULL(z, NULL); + // PS_VECTOR_CHECK_NULL(dz, NULL); psVector *zFit = NULL; @@ -548,17 +560,15 @@ // XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64) -# if (0) // moved to psLib -// XXX EAM : this version uses myPoly->nX as Norder, not Nterms -psF32 Polynomial2DEval(const psPolynomial2D* myPoly, - psF32 x, - psF32 y) -{ - PS_POLY_CHECK_NULL(myPoly, NAN); +psF64 Polynomial2DEval(const psPolynomial2D* myPoly, + psF64 x, + psF64 y) +{ + // PS_POLY_CHECK_NULL(myPoly, NAN); psS32 loop_x = 0; psS32 loop_y = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 ySum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; // XXX EAM : nX is order, not nTerms @@ -577,4 +587,2 @@ return(polySum); } -# endif - Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 4953) +++ /trunk/psphot/src/psphot.h (revision 4954) @@ -1,41 +1,35 @@ +# include +# include // for strcasecmp +# include // for unlink # include -// # include -// my additions and modifications: # include "pmObjects_EAM.h" -# include -# include -# include -# include -# include - -// # define M_PI 3.14159265358979323846264338328 /* pi */ typedef struct { - psImage *image; - psImage *mask; - psImage *noise; - psMetadata *header; + psImage *image; + psImage *mask; + psImage *noise; + psMetadata *header; } psImageData; // data to test a given PSF model type typedef struct { - pmModelType modelType; - pmPSF *psf; - psArray *sources; // pointers to the original sources - psArray *modelFLT; // model fits, floating parameters - psArray *modelPSF; // model fits, PSF parameters - psVector *mask; - psVector *metric; - psVector *fitMag; - float ApResid; - float dApResid; - float skyBias; + pmModelType modelType; + pmPSF *psf; + psArray *sources; // pointers to the original sources + psArray *modelFLT; // model fits, floating parameters + psArray *modelPSF; // model fits, PSF parameters + psVector *mask; + psVector *metric; + psVector *fitMag; + float ApResid; + float dApResid; + float skyBias; } pmPSF_Test; // structure to carry a dynamic string typedef struct { - int NLINE; - int Nline; - char *line; + int NLINE; + int Nline; + char *line; } psLine; @@ -43,10 +37,6 @@ PSPHOT_MASK_INVALID = 0x01, PSPHOT_MASK_SATURATED = 0x02, - PSPHOT_MASK_KEEP = 0x08, + PSPHOT_MASK_MARKED = 0x08, } MaskValues; - -// used by mask operations -# define AND true -# define OR false # define psMemCopy(A)(psMemIncrRefCounter((A))) @@ -116,15 +106,10 @@ psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); -//void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -//void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -//void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); -//void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); psVector *psGetRowVectorFromImage(psImage *image, psU32 row); +int psImageCountPixelMask (psImage *mask, psU8 value); // basic image functions bool psImageInit (psImage *image,...); -//void psImageSmooth (psImage *image, float sigma, float Nsigma); -//psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in); -psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius); +void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma); int psphotSaveImage (psMetadata *header, psImage *image, char *filename); @@ -135,5 +120,5 @@ // polynomial functions -psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y); +psF64 Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y); psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm); psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); @@ -143,14 +128,27 @@ psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y); + psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm); void psPolynomial1DDump (psPolynomial1D *poly); -psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly); +psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly); psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x); psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type); psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr); -// fitsource -bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y); -psMetadata *fs_args (int *argc, char **argv); -int fs_usage (); -bool DumpImage (psImage *image, char *filename); +// XXX EAM : bug : missing from psLib/*.h +bool psMinimizeGaussNewtonDelta (psVector *delta, + const psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func); + +// not included in the .h file? +bool p_psVectorPrint ( + int fd, ///< output file descriptor + psVector *a, ///< vector to print + char *name ///< name of vector (for title) +); + +bool p_psVectorPrintRow (int fd, psVector *a, char *name); Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 4953) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 4954) @@ -66,7 +66,7 @@ // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, model->radius, "or", PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); status = pmSourceFitModel (source, model, true); - psImageKeepCircle (source->mask, x, y, model->radius, "and", ~PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); if (!status || (model->params->data.F32[1] < 0)) { psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 4953) +++ /trunk/psphot/src/psphotArguments.c (revision 4954) @@ -43,5 +43,5 @@ psMetadata *config = psMetadataAlloc (); psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_META_MULTI, "folder for psf model entries", NULL); - config = psMetadataParseConfig (config, &Nfail, argv[3], FALSE); + config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE); fprintf (stderr, "loaded config...\n"); Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 4953) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 4954) @@ -50,5 +50,7 @@ // recalculate the source moments using the larger galaxy moments radius - status = pmSourceMoments (source, GAL_MOMENTS_RAD); + // XXX EAM : 0.5.0 code used _EAM version + // status = pmSourceMoments (source, GAL_MOMENTS_RAD); + status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD); if (!status) { source->type = PS_SOURCE_DROP_GALAXY; // better choice? @@ -74,7 +76,9 @@ // fit FLT (not PSF) model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + // XXX EAM : 0.5.0 code used _EAM version + // status = pmSourceFitModel (source, model, false); status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); if (!status) { // if the fit fails, we need to change the classification Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 4953) +++ /trunk/psphot/src/psphotOutput.c (revision 4954) @@ -372,18 +372,18 @@ // write the peaks to an output file -bool psPeaksWriteText (psArray *sources, char *filename) { +bool pmPeaksWriteText (psArray *sources, char *filename) { int i; FILE *f; - psPeak *peak; + pmPeak *peak; f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("psPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - peak = (psPeak *) sources->data[i]; + psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + peak = (pmPeak *) sources->data[i]; if (peak == NULL) continue; fprintf (f, "%5d %5d %7.1f\n", @@ -638,5 +638,5 @@ unlink (filename); psFits *fits = psFitsAlloc (filename); - psFitsWriteImage (fits, header, image, 0, NULL); + psFitsWriteImage (fits, header, image, 0); psFree (fits); return (TRUE); Index: /trunk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 4953) +++ /trunk/psphot/src/psphotSetup.c (revision 4954) @@ -75,5 +75,5 @@ psRegion keep = psRegionSet (XMIN, XMAX, YMIN, YMAX); keep = psRegionForImage (image, &keep); - psImageKeepRegion (mask, &keep, OR, PSPHOT_MASK_INVALID); + psImageKeepRegion (mask, &keep, "OR", PSPHOT_MASK_INVALID); // XXX EAM : psRegionForImage should not take psRegion *keep Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 4953) +++ /trunk/psphot/src/psphotSourceStats.c (revision 4954) @@ -35,5 +35,6 @@ // measure basic source moments - status = pmSourceMoments (source, RADIUS); + // XXX EAM : choose between these two versions + status = pmSourceMoments_EAM (source, RADIUS); if (!status) { psFree (source); Index: /trunk/psphot/src/pspsf.c =================================================================== --- /trunk/psphot/src/pspsf.c (revision 4953) +++ /trunk/psphot/src/pspsf.c (revision 4954) @@ -31,5 +31,5 @@ } - p_psMemSetDeallocator(psf, (psFreeFcn) pmPSFFree); + psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); return(psf); } @@ -74,5 +74,5 @@ } - p_psMemSetDeallocator(test, (psFreeFcn) pmPSF_TestFree); + psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree); return (test); } @@ -106,7 +106,7 @@ // set temporary object mask and fit object // fit model as FLT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80); + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); // exclude the poor fits @@ -142,5 +142,5 @@ y = source->peak->y; - psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80); + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); status = pmSourceFitModel (source, modelPSF, true); @@ -168,5 +168,5 @@ next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); }