Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5717) +++ /trunk/psphot/Makefile (revision 5718) @@ -39,4 +39,7 @@ $(SRC)/psphotMagnitudes.$(ARCH).o \ $(SRC)/psphotImageBackground.$(ARCH).o \ +$(SRC)/psphotBasicDeblend.$(ARCH).o \ +$(SRC)/psphotModelGroupInit.$(ARCH).o \ +$(SRC)/pmSourceContour.$(ARCH).o \ $(SRC)/psLine.$(ARCH).o \ $(SRC)/psModulesUtils.$(ARCH).o \ Index: /trunk/psphot/src/models/pmModel_TAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_TAUSS.c (revision 5718) +++ /trunk/psphot/src/models/pmModel_TAUSS.c (revision 5718) @@ -0,0 +1,164 @@ + +/****************************************************************************** + 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_TAUSS(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_TAUSS(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_TAUSS (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_TAUSS (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_TAUSS (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_TAUSS (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]; + // XXX: be wary of a bug here. EAM had his own version of psPolynomial2DEval(). + // I think the reason was that hewas using nX=xOrder. Not sure. I changed this + // without verifying. + // out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]); + out[i] = psPolynomial2DEval(poly, out[2], out[3]); + } + return(true); +} + +// check the status of the fitted model +bool pmModelFitStatus_TAUSS (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_WAUSS.c =================================================================== --- /trunk/psphot/src/models/pmModel_WAUSS.c (revision 5717) +++ /trunk/psphot/src/models/pmModel_WAUSS.c (revision 5718) @@ -59,5 +59,4 @@ 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]); @@ -66,5 +65,5 @@ } -bool psModelGuess_WAUSS (psModel *model, psSource *source) { +bool psModelGuess_WAUSS (psModel *model, pmSource *source) { psVector *params = model->params; Index: unk/psphot/src/pmModelGroup.c =================================================================== --- /trunk/psphot/src/pmModelGroup.c (revision 5717) +++ (revision ) @@ -1,111 +1,0 @@ -# include "pmModelGroup.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" -# include "models/pmModel_TGAUSS.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}, - {"PS_MODEL_TGAUSS", 8, pmModelFunc_TGAUSS, pmModelFlux_TGAUSS, pmModelRadius_TGAUSS, pmModelLimits_TGAUSS, pmModelGuess_TGAUSS, pmModelFromPSF_TGAUSS, pmModelFitStatus_TGAUSS}, -}; - -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: unk/psphot/src/pmModelGroup.h =================================================================== --- /trunk/psphot/src/pmModelGroup.h (revision 5717) +++ (revision ) @@ -1,104 +1,0 @@ - -# ifndef PM_MODEL_GROUP_H -# define PM_MODEL_GROUP_H - -# include "pmObjects_EAM.h" -# include "pmPSF.h" - -/** - * - * 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: unk/psphot/src/pmObjects_EAM.c =================================================================== --- /trunk/psphot/src/pmObjects_EAM.c (revision 5717) +++ (revision ) @@ -1,1811 +1,0 @@ -/** @file pmObjects.c - * - * This file will ... - * - * @author GLG, MHPCC - * @author EAM, IfA: significant modifications. - * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-26 03:00:41 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#include -#include -#include -#include "pslib.h" -#include "psConstants.h" -#include "pmObjects_EAM.h" -#include "pmModelGroup.h" - - -/****************************************************************************** -pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. -*****************************************************************************/ -pmPeak *pmPeakAlloc(psS32 x, - psS32 y, - psF32 counts, - pmPeakType class) -{ - pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); - tmp->x = x; - tmp->y = y; - tmp->counts = counts; - tmp->class = class; - - return(tmp); -} - -/****************************************************************************** -pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members -to zero. -*****************************************************************************/ -pmMoments *pmMomentsAlloc() -{ - pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); - tmp->x = 0.0; - tmp->y = 0.0; - tmp->Sx = 0.0; - tmp->Sx = 0.0; - tmp->Sxy = 0.0; - tmp->Sum = 0.0; - tmp->Peak = 0.0; - tmp->Sky = 0.0; - tmp->nPixels = 0; - - return(tmp); -} - -static void modelFree(pmModel *tmp) -{ - psFree(tmp->params); - psFree(tmp->dparams); -} - -/****************************************************************************** -pmModelAlloc(): Allocate the pmModel structure, along with its parameters, -and initialize the type member. Initialize the params to 0.0. -XXX EAM: simplifying code with pmModelParameterCount -*****************************************************************************/ -pmModel *pmModelAlloc(pmModelType type) -{ - pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); - - tmp->type = type; - tmp->chisq = 0.0; - tmp->nIter = 0; - psS32 Nparams = pmModelParameterCount (type); - if (Nparams == 0) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NULL); - } - - tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); - - for (psS32 i = 0; i < tmp->params->n; i++) { - tmp->params->data.F32[i] = 0.0; - tmp->dparams->data.F32[i] = 0.0; - } - - psMemSetDeallocator(tmp, (psFreeFunc) modelFree); - return(tmp); -} - -/****************************************************************************** -XXX EAM : we can now free these pixels - memory ref is incremented now -*****************************************************************************/ -static void sourceFree(pmSource *tmp) -{ - psFree(tmp->peak); - psFree(tmp->pixels); - psFree(tmp->weight); - psFree(tmp->mask); - psFree(tmp->moments); - psFree(tmp->modelPSF); - psFree(tmp->modelFLT); -} - -/****************************************************************************** -pmSourceAlloc(): Allocate the pmSource structure and initialize its members -to NULL. -*****************************************************************************/ -pmSource *pmSourceAlloc() -{ - pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); - tmp->peak = NULL; - tmp->pixels = NULL; - tmp->weight = NULL; - tmp->mask = NULL; - tmp->moments = NULL; - tmp->modelPSF = NULL; - tmp->modelFLT = NULL; - tmp->type = 0; - psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); - - return(tmp); -} - -/****************************************************************************** -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. - -XXX: What types should be supported? Only F32 is implemented. - -XXX: We currently step through the input vector twice; once to determine the -size of the output vector, then to set the values of the output vector. -Depending upon actual use, this may need to be optimized. -*****************************************************************************/ -psVector *pmFindVectorPeaks(const psVector *vector, - psF32 threshold) -{ - PS_ASSERT_VECTOR_NON_NULL(vector, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); - PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL); - int count = 0; - int n = vector->n; - - // - // Special case: the input vector has a single element. - // - if (n == 1) { - psVector *tmpVector = NULL; - ; - if (vector->data.F32[0] > threshold) { - tmpVector = psVectorAlloc(1, PS_TYPE_U32); - tmpVector->data.U32[0] = 0; - } else { - tmpVector = psVectorAlloc(0, PS_TYPE_U32); - } - return(tmpVector); - } - - // - // Determine if first pixel is a peak - // - if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { - count++; - } - - // - // Determine if interior pixels are peaks - // - 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)) { - count++; - } - } - - // - // Determine if last pixel is a peak - // - if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { - count++; - } - - // - // We know how many peaks exist, so we now allocate a psVector to store - // those peaks. - // - psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32); - count = 0; - - // - // Determine if first pixel is a peak - // - if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { - tmpVector->data.U32[count++] = 0; - } - - // - // Determine if interior pixels are peaks - // - 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)) { - tmpVector->data.U32[count++] = i; - } - } - - // - // Determine if last pixel is a peak - // - if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { - tmpVector->data.U32[count++] = n-1; - } - - return(tmpVector); -} - -/****************************************************************************** -getRowVectorFromImage(): a private function which simply returns a -psVector containing the specified row of data from the psImage. - -XXX: Is there a better way to do this? -XXX EAM: does this really need to alloc a new vector??? -*****************************************************************************/ -static psVector *getRowVectorFromImage(psImage *image, - psU32 row) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - for (psU32 col = 0; col < image->numCols ; col++) { - tmpVector->data.F32[col] = image->data.F32[row][col]; - } - return(tmpVector); -} - -/****************************************************************************** -myListAddPeak(): A private function which allocates a psArray, if the list -argument is NULL, otherwise it adds the peak to that list. -XXX EAM : changed the output to psArray -XXX EAM : Switched row, col args -XXX EAM : NOTE: this was changed in the call, so the new code is consistent -*****************************************************************************/ -static psArray *myListAddPeak(psArray *list, - psS32 row, - psS32 col, - psF32 counts, - pmPeakType type) -{ - pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); - - if (list == NULL) { - list = psArrayAlloc(100); - list->n = 0; - } - psArrayAdd(list, 100, tmpPeak); - psFree (tmpPeak); - // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) - - return(list); -} - -/****************************************************************************** -pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage -above the given threshold. Returns a psArray containing location (x/y value) -of all peaks. - -XXX: I'm not convinced the peak type definition in the SDRS is mutually -exclusive. Some peaks can have multiple types. Edges for sure. Also, a -digonal line with the same value at each point will have a peak for every -point on that line. - -XXX: This does not work if image has either a single row, or a single column. - -XXX: In the output psArray elements, should we use the image row/column offsets? - Currently, we do not. - -XXX: Merge with CVS 1.20. This had the proper code for images with a single -row or column. -*****************************************************************************/ -psArray *pmFindImagePeaks(const psImage *image, - psF32 threshold) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - if ((image->numRows == 1) || (image->numCols == 1)) { - psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); - return(NULL); - } - psVector *tmpRow = NULL; - psU32 col = 0; - psU32 row = 0; - psArray *list = NULL; - - // - // Find peaks in row 0 only. - // - row = 0; - tmpRow = getRowVectorFromImage((psImage *) image, row); - psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); - - for (psU32 i = 0 ; i < row1->n ; i++ ) { - col = row1->data.U32[i]; - // - // Determine if pixel (0,0) is a peak. - // - 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])) { - - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } 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])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } 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])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); - } - } - psFree (tmpRow); - psFree (row1); - - // - // Exit if this image has a single row. - // - if (image->numRows == 1) { - return(list); - } - - // - // Find peaks in interior rows only. - // - for (row = 1 ; row < (image->numRows - 1) ; row++) { - tmpRow = getRowVectorFromImage((psImage *) image, row); - row1 = pmFindVectorPeaks(tmpRow, threshold); - - // Step through all local peaks in this row. - for (psU32 i = 0 ; i < row1->n ; i++ ) { - pmPeakType myType = PM_PEAK_UNDEF; - col = row1->data.U32[i]; - - if (col == 0) { - // 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])) { - myType = PM_PEAK_EDGE; - list = myListAddPeak(list, row, col, image->data.F32[row][col], myType); - } - } else if (col < (image->numCols - 1)) { - // 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])) { - 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])) { - 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])) { - myType = PM_PEAK_FLAT; - } - - list = myListAddPeak(list, row, col, image->data.F32[row][col], myType); - } - } - } else if (col == (image->numCols - 1)) { - // 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])) { - myType = PM_PEAK_EDGE; - list = myListAddPeak(list, row, col, image->data.F32[row][col], myType); - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); - } - - } - psFree (tmpRow); - psFree (row1); - } - - // - // Find peaks in the last row only. - // - row = image->numRows - 1; - tmpRow = getRowVectorFromImage((psImage *) image, row); - row1 = pmFindVectorPeaks(tmpRow, threshold); - for (psU32 i = 0 ; i < row1->n ; i++ ) { - col = row1->data.U32[i]; - 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])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } 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])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } 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])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range."); - } - } - psFree (tmpRow); - psFree (row1); - return(list); -} - -// XXX: Macro this. -static bool isItInThisRegion(const psRegion *valid, - psS32 x, - psS32 y) -{ - - // XXX EAM this function should check the status of valid - if (valid == NULL) return (true); - - if ((x >= valid->x0) && - (x <= valid->x1) && - (y >= valid->y0) && - (y <= valid->y1)) { - return(true); - } - - return(false); -} - - -/****************************************************************************** -psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have -a peak value above the given maximum, or fall outside the valid region. - -XXX: Should the sky value be used when comparing the maximum? - -XXX: warning message if valid is NULL? - -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, - psF32 maxValue, - const psRegion *valid) -{ - PS_ASSERT_PTR_NON_NULL(peaks, NULL); - - psListElem *tmpListElem = (psListElem *) peaks->head; - psS32 indexNum = 0; - - // printf("pmCullPeaks(): list size is %d\n", peaks->size); - while (tmpListElem != NULL) { - pmPeak *tmpPeak = (pmPeak *) tmpListElem->data; - if ((tmpPeak->counts > maxValue) || - ((valid != NULL) && - (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { - psListRemoveData(peaks, (psPtr) tmpPeak); - } - - indexNum++; - tmpListElem = tmpListElem->next; - } - - return(peaks); -} - -// XXX EAM: I changed this to return a new, subset array -// rather than alter the existing one -psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid) -{ - PS_ASSERT_PTR_NON_NULL(peaks, NULL); - - psArray *output = psArrayAlloc (200); - output->n = 0; - - psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); - - for (int i = 0; i < peaks->n; i++) { - pmPeak *tmpPeak = (pmPeak *) peaks->data[i]; - if (tmpPeak->counts > maxValue) - continue; - if (valid != NULL) { - if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) - continue; - } - psArrayAdd (output, 200, tmpPeak); - } - return(output); -} - -/****************************************************************************** -pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this -routine creates a new pmSource data structure and sets the following members: - ->pmPeak - ->pmMoments->sky - -The sky value is set from the pixels in the square annulus surrounding the -peak pixel. - -We simply create a subSet image and mask the inner pixels, then call -psImageStats on that subImage+mask. - -XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. - -XXX: Use static data structures for: - subImage - subImageMask - myStats - -XXX: ensure that the inner and out radius fit in the actual image. Should - we generate an error, or warning? Currently an error. - -XXX: Sync with IfA on whether the peak x/y coords are data structure coords, - or they use the image row/column offsets. - -XXX: Should we simply set pmSource->peak = peak? If so, should we increase -the reference counter? Or, should we copy the data structure? - -XXX: Currently the subimage always has an even number of rows/columns. Is - this correct? Since there is a center pixel, maybe it should have an - odd number of rows/columns. - -XXX: Use psTrace() for the print statements. - -XXX: Don't use separate structs for the subimage and mask. Use the source-> - members. -*****************************************************************************/ -bool pmSourceLocalSky (pmSource *source, - psStatsOptions statsOptions, - psF32 Radius) -{ - - psImage *image = source->pixels; - psImage *mask = source->mask; - pmPeak *peak = source->peak; - psRegion srcRegion; - - // XXX EAM : psRegionXXX funcs should take 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", ~PSPHOT_MASK_MARKED); - - psF64 tmpF64; - p_psGetStatValue(myStats, &tmpF64); - psFree (myStats); - - if (isnan(tmpF64)) return (false); - source->moments = pmMomentsAlloc(); - source->moments->Sky = (psF32) tmpF64; - return (true); -} - -/****************************************************************************** -bool checkRadius2(): private function which simply determines if the (x, y) -point is within the radius of the specified peak. - -XXX: macro this for performance. -XXX: this is rather inefficient - at least compute and compare against radius^2 -*****************************************************************************/ -static bool checkRadius2(psF32 xCenter, - psF32 yCenter, - psF32 radius, - psF32 x, - psF32 y) -{ - /// XXX EAM should compare with hypot (x,y) for speed - if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { - return(true); - } - - return(false); -} - -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -pmSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - pmSource - pmSource->peak - pmSource->pixels - pmSource->weight - pmSource->mask - -XXX: The peak calculations are done in image coords, not subImage coords. - -XXX EAM : this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons -*****************************************************************************/ -# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -bool pmSourceMoments(pmSource *source, - psF32 radius) -{ - // 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->weight, NULL); - PS_FLOAT_COMPARE(0.0, radius, NULL); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - psF32 xPeak = source->peak->x; - psF32 yPeak = source->peak->y; - - // XXX why do I get different results for these two methods of finding Sx? - // XXX Sx, Sy would be better measured if we clip pixels close to sky - // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? - // We loop through all pixels in this subimage (source->pixels), and for each - // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col])) { - continue; - } - - psF32 xDiff = col + source->pixels->col0 - xPeak; - psF32 yDiff = row + source->pixels->row0 - yPeak; - - // XXX EAM : calculate xDiff, yDiff up front; - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) { - continue; - } - - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - - // XXX EAM : check for valid S/N in pixel - // XXX EAM : should this limit be user-defined? - if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) { - continue; - } - - Sum += pDiff; - X1 += xDiff * pDiff; - Y1 += yDiff * pDiff; - XY += xDiff * yDiff * pDiff; - - X2 += PS_SQR(xDiff) * pDiff; - Y2 += PS_SQR(yDiff) * pDiff; - - peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); - numPixels++; - } - } - // 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, - "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace (".psModules.pmSourceMoments", 5, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - return (false); - } - - source->moments->x = x + xPeak; - source->moments->y = y + yPeak; - - // XXX EAM : Sxy needs to have x*y subtracted - source->moments->Sxy = XY/Sum - x*y; - source->moments->Sum = Sum; - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - - // XXX EAM : these values can be negative, so we need to limit the range - source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); - source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); - - psTrace (".psModules.pmSourceMoments", 4, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, source->moments->Sxy); - - return(true); -} - -// XXX EAM : I used -int pmComparePeakAscend (const void **a, const void **b) -{ - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) - return (-1); - if (diff > FLT_EPSILON) - return (+1); - return (0); -} - -int pmComparePeakDescend (const void **a, const void **b) -{ - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) - return (+1); - if (diff > FLT_EPSILON) - return (-1); - return (0); -} - -/****************************************************************************** -pmSourcePSFClump(source, metadata): Find the likely PSF clump in the -sigma-x, sigma-y plane. return 0,0 clump in case of error. -*****************************************************************************/ - -// XXX EAM include a S/N cutoff in selecting the sources? -pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) -{ - -# define NPIX 10 -# define SCALE 0.1 - - psArray *peaks = NULL; - pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; - pmPSFClump psfClump = emptyClump; - - PS_ASSERT_PTR_NON_NULL(sources, emptyClump); - PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); - - // find the sigmaX, sigmaY clump - { - psStats *stats = NULL; - psImage *splane = NULL; - int binX, binY; - - // construct a sigma-plane image - // psImageAlloc does zero the data - 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)); - } - - // place the sources in the sigma-plane image (ignore 0,0 values?) - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - if (tmpSrc == NULL) { - continue; - } - if (tmpSrc->moments == NULL) { - continue; - } - - // Sx,Sy are limited at 0. a peak at 0,0 is artificial - if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { - continue; - } - - // for the moment, force splane dimensions to be 10x10 image pix - binX = tmpSrc->moments->Sx/SCALE; - if (binX < 0) - continue; - if (binX >= splane->numCols) - continue; - - binY = tmpSrc->moments->Sy/SCALE; - if (binY < 0) - continue; - if (binY >= splane->numRows) - continue; - - splane->data.F32[binY][binX] += 1.0; - } - - // find the peak in this image - stats = psStatsAlloc (PS_STAT_MAX); - stats = psImageStats (stats, splane, NULL, 0); - peaks = pmFindImagePeaks (splane, stats[0].max / 2); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - - psFree (splane); - psFree (stats); - - } - // XXX EAM : possible errors: - // 1) no peak in splane - // 2) no significant peak in splane - - // measure statistics on Sx, Sy if Sx, Sy within range of clump - { - pmPeak *clump; - psF32 minSx, maxSx; - psF32 minSy, maxSy; - psVector *tmpSx = NULL; - psVector *tmpSy = NULL; - psStats *stats = NULL; - - // XXX EAM : this lets us takes the single highest peak - psArraySort (peaks, pmComparePeakDescend); - clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); - - // define section window for clump - minSx = clump->x * SCALE - 0.2; - maxSx = clump->x * SCALE + 0.2; - minSy = clump->y * SCALE - 0.2; - maxSy = clump->y * SCALE + 0.2; - - tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSx->n = 0; - tmpSy->n = 0; - - // XXX clip sources based on flux? - // create vectors with Sx, Sy values in window - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - if (tmpSrc->moments->Sx < minSx) - continue; - if (tmpSrc->moments->Sx > maxSx) - continue; - if (tmpSrc->moments->Sy < minSy) - continue; - if (tmpSrc->moments->Sy > maxSy) - continue; - tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; - tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; - tmpSx->n++; - tmpSy->n++; - if (tmpSx->n == tmpSx->nalloc) { - psVectorRealloc (tmpSx, tmpSx->nalloc + 100); - psVectorRealloc (tmpSy, tmpSy->nalloc + 100); - } - } - - // measures stats of Sx, Sy - stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - - stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - psfClump.X = stats->clippedMean; - psfClump.dX = stats->clippedStdev; - - stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); - psfClump.Y = stats->clippedMean; - psfClump.dY = stats->clippedStdev; - - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); - // these values should be pushed on the metadata somewhere - - psFree (stats); - psFree (peaks); - psFree (tmpSx); - psFree (tmpSy); - } - - return (psfClump); -} - -/****************************************************************************** -pmSourceRoughClass(source, metadata): make a guess at the source -classification. - -XXX: push the clump info into the metadata? - -XXX: How can this function ever return FALSE? - -XXX EAM : add the saturated mask value to metadata -*****************************************************************************/ - -bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) -{ - - psBool rc = true; - - int Nsat = 0; - int Ngal = 0; - int Nstar = 0; - int Npsf = 0; - int Ncr = 0; - int Nsatstar = 0; - psRegion allArray = psRegionSet (0, 0, 0, 0); - - psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); - starsn->n = 0; - - // check return status value (do these exist?) - bool status; - psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); - psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); - psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); - // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); - - // XXX allow clump size to be scaled relative to sigmas? - // make rough IDs based on clumpX,Y,DX,DY - for (psS32 i = 0 ; i < sources->n ; i++) { - - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - tmpSrc->peak->class = 0; - - psF32 sigX = tmpSrc->moments->Sx; - psF32 sigY = tmpSrc->moments->Sy; - - // calculate and save signal-to-noise estimates - psF32 S = tmpSrc->moments->Sum; - psF32 A = 4 * M_PI * sigX * sigY; - psF32 B = tmpSrc->moments->Sky; - 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 SATURATE if mask is NULL? - int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED); - - // saturated star (size consistent with PSF or larger) - // Nsigma should be user-configured parameter - bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); - if ((Nsatpix > 1) && big) { - tmpSrc->type |= PM_SOURCE_SATSTAR; - Nsatstar ++; - continue; - } - - // saturated object (not a star, eg bleed trails, hot pixels) - if (Nsatpix > 1) { - tmpSrc->type |= PM_SOURCE_SATURATED; - Nsat ++; - continue; - } - - // likely defect (too small to be stellar) (push out to 3 sigma) - // low S/N objects which are small are probably stellar - // only set candidate defects if - if ((sigX < 0.05) || (sigY < 0.05)) { - tmpSrc->type |= PM_SOURCE_DEFECT; - Ncr ++; - continue; - } - - // likely unsaturated galaxy (too large to be stellar) - if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - tmpSrc->type |= PM_SOURCE_GALAXY; - Ngal ++; - continue; - } - - // the rest are probable stellar objects - starsn->data.F32[starsn->n] = SN; - starsn->n ++; - Nstar ++; - - // PSF star (within 1.5 sigma of clump center - psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); - if ((SN > PSF_SN_LIM) && (radius < 1.5)) { - tmpSrc->type |= PM_SOURCE_PSFSTAR; - Npsf ++; - continue; - } - - // random type of star - tmpSrc->type |= PM_SOURCE_OTHER; - } - - { - psStats *stats = NULL; - stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - stats = psVectorStats (stats, starsn, NULL, NULL, 0); - psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); - psFree (stats); - psFree (starsn); - } - - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); - - return(rc); -} - -/****************************************************************************** -pmSourceSetPixelsCircle(source, image, radius) - -XXX: Why boolean output? - -XXX: Why are we checking source->moments for NULL? Should the circle be - centered on the centroid or the peak? - -XXX: The circle will have a diameter of (1+radius). This is different from - the pmSourceSetLocal() function. - -XXX: this should probably use the psImageMaskCircle Function for consistency -*****************************************************************************/ -bool pmSourceSetPixelsCircle(pmSource *source, - const psImage *image, - psF32 radius) -{ - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_FLOAT_COMPARE(0.0, radius, false); - - // - // We define variables for code readability. - // - // XXX: Since the peak->xy coords are in image, not subImage coords, - // these variables should be renamed for clarity (imageCenterRow, etc). - // - psS32 radiusS32 = (psS32) radius; - psS32 SubImageCenterRow = source->peak->y; - psS32 SubImageCenterCol = source->peak->x; - // XXX EAM : for the circle to stay on the image - // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) - psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); - psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); - psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); - psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); - - // XXX: Must recycle image. - // XXX EAM: this message reflects a programming error we know about. - // i am setting it to a trace message which we can take out - if (source->pixels != NULL) { - psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, - "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); - psFree(source->pixels); - } - source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, - SubImageStartRow, - SubImageEndCol, - SubImageEndRow)); - - // XXX: Must recycle image. - if (source->mask != NULL) { - psFree(source->mask); - } - source->mask = psImageAlloc(source->pixels->numCols, - source->pixels->numRows, - PS_TYPE_U8); // XXX EAM : type was F32 - - // - // Loop through the subimage mask, initialize mask to 0 or 1. - // XXX EAM: valid pixels should have 0, not 1 - for (psS32 row = 0 ; row < source->mask->numRows; row++) { - for (psS32 col = 0 ; col < source->mask->numCols; col++) { - - if (checkRadius2((psF32) radiusS32, - (psF32) radiusS32, - radius, - (psF32) col, - (psF32) row)) { - source->mask->data.U8[row][col] = 0; - } else { - source->mask->data.U8[row][col] = 1; - } - } - } - return(true); -} - -/****************************************************************************** -pmSourceModelGuess(source, model): This function allocates a new -pmModel structure based on the given modelType specified in the argument list. -The corresponding pmModelGuess function is returned, and used to -supply the values of the params array in the pmModel structure. - -XXX: Many of the initial parameters are set to 0.0 since I don't know what -the appropiate initial guesses are. - -XXX: Many parameters are based on the src->moments structure, which is in -image, not subImage coords. Therefore, the calls to the model evaluation -functions will be in image, not subImage coords. Remember this. -*****************************************************************************/ -pmModel *pmSourceModelGuess(pmSource *source, - pmModelType modelType) -{ - 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); -} - -/****************************************************************************** -evalModel(source, level, row): a private function which evaluates the -source->modelPSF function at the specified coords. The coords are subImage, not -image coords. - -NOTE: The coords are in subImage source->pixel coords, not image coords. - -XXX: reverse order of row,col args? - -XXX: rename all coords in this file such that their name defines whether -the coords is in subImage or image space. - -XXX: This should probably be a public pmModules function. - -XXX: Use static vectors for x. - -XXX: Figure out if it's (row, col) or (col, row) for the model functions. - -XXX: For a while, the first psVectorAlloc() was generating a seg fault during -testing. Try to reproduce that and debug. -*****************************************************************************/ - -// 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. - // - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - x->data.F32[0] = (psF32) (col + image->col0); - x->data.F32[1] = (psF32) (row + image->row0); - psF32 tmpF; - pmModelFunc modelFunc; - - modelFunc = pmModelFunc_GetFunction (model->type); - tmpF = modelFunc (NULL, model->params, x); - psFree(x); - return(tmpF); -} - -/****************************************************************************** -findValue(source, level, row, col, dir): a private function which determines -the column coordinate of the model function which has the value "level". If -dir equals 0, then you loop leftwards from the peak pixel, otherwise, -rightwards. - -XXX: reverse order of row,col args? - -XXX: Input row/col are in image coords. - -XXX: The result is returned in image coords. -*****************************************************************************/ -static psF32 findValue(pmSource *source, - psF32 level, - psU32 row, - psU32 col, - psU32 dir) -{ - // - // Convert coords to subImage space. - // - psU32 subRow = row - source->pixels->row0; - psU32 subCol = col - source->pixels->col0; - - // Ensure that the starting column is allowable. - if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { - psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); - return(NAN); - } - if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { - psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); - return(NAN); - } - - // XXX EAM : i changed this to match pmModelEval above, but see - // XXX EAM the note below in pmSourceContour - psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); - if (oldValue == level) { - return(((psF32) (subCol + source->pixels->col0))); - } - - // - // We define variables incr and lastColumn so that we can use the same loop - // whether we are stepping leftwards, or rightwards. - // - psS32 incr; - psS32 lastColumn; - if (dir == 0) { - incr = -1; - lastColumn = -1; - } else { - incr = 1; - lastColumn = source->pixels->numCols; - } - subCol+=incr; - - while (subCol != lastColumn) { - psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); - if (oldValue == level) { - return((psF32) (subCol + source->pixels->col0)); - } - - if ((newValue <= level) && (level <= oldValue)) { - // This is simple linear interpolation. - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); - } - - if ((oldValue <= level) && (level <= newValue)) { - // This is simple linear interpolation. - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); - } - - subCol+=incr; - } - - return(NAN); -} - -/****************************************************************************** -pmSourceContour(src, img, level, mode): For an input subImage, and model, this -routine returns a psArray of coordinates that evaluate to the specified level. - -XXX: Probably should remove the "image" argument. -XXX: What type should the output coordinate vectors consist of? col,row? -XXX: Why a pmArray output? -XXX: doex x,y correspond with col,row or row/col? -XXX: What is mode? -XXX: The top, bottom of the contour is not correctly determined. -XXX EAM : this functions is using the model for the contour, but it should - be using only the image counts -*****************************************************************************/ -psArray *pmSourceContour(pmSource *source, - const psImage *image, - psF32 level, - pmContourType mode) -{ - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(image, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->modelFLT, false); - // XXX EAM : what is the purpose of modelPSF/modelFLT? - - // - // Allocate data for x/y pairs. - // - psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - - // - // Start at the row with peak pixel, then decrement. - // - psS32 col = source->peak->x; - for (psS32 row = source->peak->y; row>= 0 ; row--) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - //printf("HERE the intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - // - // Start at the row (+1) with peak pixel, then increment. - // - col = source->peak->x; - for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - - // - // Allocate an array for result, store coord vectors there. - // - psArray *tmpArray = psArrayAlloc(2); - tmpArray->data[0] = (psPtr *) yVec; - tmpArray->data[1] = (psPtr *) xVec; - return(tmpArray); -} - -// XXX EAM : these are better starting values, but should be available from metadata? -#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 -#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 -/****************************************************************************** -pmSourceFitModel(source, model): must create the appropiate arguments to the -LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. - -XXX: should there be a mask value? -XXX EAM : fit the specified model (not necessarily the one in source) -*****************************************************************************/ -bool pmSourceFitModel_v5(pmSource *source, - pmModel *model, - const bool PSF) -{ - 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->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - - // XXX EAM : is it necessary for the mask & weight to exist? the - // tests below could be conditions (!NULL) - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - int nParams = PSF ? params->n - 4 : params->n; - - // find the number of valid pixels - // XXX EAM : this loop and the loop below could just be one pass - // using the psArrayAdd and psVectorExtend functions - psS32 count = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - count++; - } - } - } - if (count < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - return(false); - } - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(count); - psVector *y = psVectorAlloc(count, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); - psS32 tmpCnt = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - // XXX: Convert i/j to image space: - // XXX EAM: coord order is (x,y) == (col,row) - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[tmpCnt] = (psPtr *) coord; - y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; - yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]); - // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then - // the minimization function calculates sq() - tmpCnt++; - } - } - } - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - 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; - } - } - - // XXX EAM : covar must be F64? - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); - - psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); - fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if ((paramMask != NULL) && paramMask->data.U8[i]) continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - - // XXX EAM: we need to do something (give an error?) if rc is false - // XXX EAM: psMinimizeLMChi2 does not check convergence - - // XXX models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - - // XXX EAM: save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // XXX EAM get the Gauss-Newton distance for fixed model parameters - if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if (!paramMask->data.U8[i]) continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - psFree (delta); - } - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myMin); - psFree(covar); - psFree(paramMask); - - rc = (onPic && fitStatus); - return(rc); -} - -// XXX EAM : new version with parameter range limits and weight enhancement -bool pmSourceFitModel (pmSource *source, - pmModel *model, - const bool PSF) -{ - 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->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - // XXX EAM : is it necessary for the mask & weight to exist? the - // tests below could be conditions (!NULL) - - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - psF32 Ro, ymodel; - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - // XXX EAM : I need to use the sky value to constrain the weight model - int nParams = PSF ? params->n - 4 : params->n; - psF32 So = params->data.F32[0]; - - // find the number of valid pixels - // XXX EAM : this loop and the loop below could just be one pass - // using the psArrayAdd and psVectorExtend functions - psS32 count = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - count++; - } - } - } - if (count < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - return(false); - } - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(count); - psVector *y = psVectorAlloc(count, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); - psS32 tmpCnt = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - // XXX: Convert i/j to image space: - // XXX EAM: coord order is (x,y) == (col,row) - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[tmpCnt] = (psPtr *) coord; - y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; - - // compare observed flux to model flux to adjust weight - ymodel = modelFunc (NULL, model->params, coord); - - // this test modifies the weight based on deviation from the model flux - Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); - - // psMinimizeLMChi2 takes wt = 1/dY^2 - if (source->weight->data.F32[i][j] == 0) { - yErr->data.F32[tmpCnt] = 0.0; - } else { - yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro); - } - tmpCnt++; - } - } - } - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - 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; - } - } - - - // XXX EAM : I've added three types of parameter range checks - // XXX EAM : this requires my new psMinimization functions - pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); - psVector *beta_lim = NULL; - psVector *params_min = NULL; - psVector *params_max = NULL; - - // XXX EAM : in this implementation, I pass in the limits with the covar matrix. - // in the SDRS, I define a new psMinimization which will take these in - psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); - modelLimits (&beta_lim, ¶ms_min, ¶ms_max); - for (int i = 0; i < params->n; i++) { - covar->data.F64[0][i] = beta_lim->data.F32[i]; - covar->data.F64[1][i] = params_min->data.F32[i]; - covar->data.F64[2][i] = params_max->data.F32[i]; - } - - psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); - - fitStatus = psMinimizeLMChi2 (myMin, covar, params, paramMask, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if ((paramMask != NULL) && paramMask->data.U8[i]) continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - // XXX EAM: we need to do something (give an error?) if rc is false - // XXX EAM: psMinimizeLMChi2 does not check convergence - - // XXX models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - - // XXX EAM: save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // XXX EAM get the Gauss-Newton distance for fixed model parameters - if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if (!paramMask->data.U8[i]) continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - } - - psFree(paramMask); - psFree(x); - psFree(y); - psFree(myMin); - - rc = (onPic && fitStatus); - return(rc); -} - -bool p_pmSourceAddOrSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - psS32 flag) -{ - - 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; - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psS32 imageCol; - psS32 imageRow; - - for (psS32 i = 0; i < image->numRows; i++) { - for (psS32 j = 0; j < image->numCols; j++) { - if ((mask != NULL) && mask->data.U8[i][j]) continue; - psF32 pixelValue; - // XXX: Should you be adding the pixels for the entire subImage, - // or a radius of pixels around it? - - // Convert i/j to imace coord space: - // XXX: Make sure you have col/row order correct. - // XXX EAM : 'center' option changes this - // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] - if (center) { - imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; - imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; - } else { - imageCol = j + image->col0; - imageRow = i + image->row0; - } - - x->data.F32[0] = (float) imageCol; - x->data.F32[1] = (float) imageRow; - pixelValue = modelFunc (NULL, params, x); - // fprintf (stderr, "%f %f %d %d %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue); - - if (flag == 1) { - pixelValue = -pixelValue; - } - - // XXX: Must figure out how to calculate the image coordinates and - // how to use the boolean "center" flag. - - image->data.F32[i][j]+= pixelValue; - } - } - psFree(x); - return(true); -} - - - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceAddModel(psImage *image, - psImage *mask, - pmModel *model, - bool center) -{ - return(p_pmSourceAddOrSubModel(image, mask, model, center, 0)); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center) -{ - return(p_pmSourceAddOrSubModel(image, mask, model, center, 1)); -} Index: unk/psphot/src/pmObjects_EAM.h =================================================================== --- /trunk/psphot/src/pmObjects_EAM.h (revision 5717) +++ (revision ) @@ -1,313 +1,0 @@ -/** @file pmObjects.h - * - * This file will ... - * - * @author GLG, MHPCC - * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-26 03:00:41 $ - * - * 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" -#include "psEllipse.h" - -# define psMemCopy(A)(psMemIncrRefCounter((A))) - -// defined mask values for mask image pixels. -// this is probably a bad solution: we will want to set mask values -// outside of the PSPHOT code. Perhaps we can set up a registered -// set of mask values with specific meanings that other functions -// can add to or define? -enum { - PSPHOT_MASK_CLEAR = 0x00, - PSPHOT_MASK_INVALID = 0x01, - PSPHOT_MASK_SATURATED = 0x02, - PSPHOT_MASK_MARKED = 0x08, -} psphotMaskValues; - -/** 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; - -/** pmPSFClump data structure - * - * - * - */ -typedef struct -{ - float X; - float dX; - float Y; - float dY; -} pmPSFClump; - -typedef int pmModelType; - -/** 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 { - PM_SOURCE_DEFECT, - PM_SOURCE_SATURATED, - - PM_SOURCE_SATSTAR, - PM_SOURCE_PSFSTAR, - PM_SOURCE_GOODSTAR, - - PM_SOURCE_POOR_FIT_PSF, - PM_SOURCE_FAIL_FIT_PSF, - PM_SOURCE_FAINTSTAR, - - PM_SOURCE_GALAXY, - PM_SOURCE_FAINT_GALAXY, - PM_SOURCE_DROP_GALAXY, - PM_SOURCE_FAIL_FIT_GAL, - PM_SOURCE_POOR_FIT_GAL, - - PM_SOURCE_OTHER, -} pmSourceType; - -/** 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 *weight; ///< Image variance. - 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. - float apMag; - float fitMag; -} -pmSource; - -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): - -*****************************************************************************/ -bool pmSourceLocalSky( - pmSource *source, ///< The input image (float) - psStatsOptions statsOptions, ///< The statistic used in calculating the background sky - float Radius ///< The inner radius of the square annulus to exclude - ); - -/****************************************************************************** - *****************************************************************************/ -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 pmSource - psMetadata *metadata ///< Contains classification parameters -); - -/****************************************************************************** -pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the -source classification. -*****************************************************************************/ -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 -); - -bool pmSourceFitModel_v5( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or FLT? -); - -bool pmSourceFitModel_v7( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or FLT? -); - -#endif Index: unk/psphot/src/pmPSF.c =================================================================== --- /trunk/psphot/src/pmPSF.c (revision 5717) +++ (revision ) @@ -1,127 +1,0 @@ -# include -# include "psLibUtils.h" -# include "pmObjects_EAM.h" -# include "pmPSF.h" -# include "pmModelGroup.h" - -/******** pmPSF functions **************************************************/ - -// function to free a pmPSF structure -static void pmPSFFree (pmPSF *psf) { - - if (psf == NULL) return; - - psFree (psf->params); - return; -} - -// allocate a pmPSF -// a PSF always has 4 parameters fewer than the equivalent model -pmPSF *pmPSFAlloc (pmModelType type) { - - int Nparams; - - pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF)); - - psf->type = type; - psf->chisq = 0.0; - psf->ApResid = 0.0; - psf->dApResid = 0.0; - psf->skyBias = 0.0; - - Nparams = pmModelParameterCount (type); - if (!Nparams) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NULL); - } - - psf->params = psArrayAlloc (Nparams - 4); - for (int i = 0; i < psf->params->n; i++) { - // XXX EAM : make polynomial order user-specified? - if (1) { - psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); - } else { - psf->params->data[i] = psPolynomial2DAlloc(0, 0, PS_POLYNOMIAL_ORD); - } - } - - psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); - return(psf); -} - -// build a PSF from a collection of free-fitted models -// the PSF ignores the first 4 (independent) model parameters -// and constructs a polynomial fit to the remaining as a function -// of image coordinate -// input: an array of pmModels, pre-allocated psf -// some of the array entries may be NULL (failed fits), ignore them -bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) { - - // construct the fit vectors from the collection of objects - // use the mask to ignore missing fits - psVector *x = psVectorAlloc (models->n, PS_TYPE_F64); - psVector *y = psVectorAlloc (models->n, PS_TYPE_F64); - psVector *z = psVectorAlloc (models->n, PS_TYPE_F64); - psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64); - - for (int i = 0; i < models->n; i++) { - pmModel *model = models->data[i]; - if (model == NULL) continue; - - // XXX EAM : this is fragile: x and y must be stored consistently at 2,3 - // note that the data is provided in the F64 array - x->data.F64[i] = model->params->data.F32[2]; - y->data.F64[i] = model->params->data.F32[3]; - } - - // we are doing a robust fit. after each pass, we drop points which are - // more deviant than three sigma. - // mask is currently updated for each pass. this is inconsistent? - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - - for (int i = 0; i < psf->params->n; i++) { - for (int j = 0; j < models->n; j++) { - pmModel *model = models->data[j]; - if (model == NULL) continue; - z->data.F64[j] = model->params->data.F32[i + 4]; - dz->data.F64[j] = 1; - // XXX EAM : need to use actual errors? - // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4 - } - - // XXX EAM : this is the expected API (cycle 7? cycle 8?) - psf->params->data[i] = psVectorClipFitPolynomial2D (psf->params->data[i], stats, mask, 0xff, z, dz, x, y); - - // XXX EAM : drop this when above is implemented... - // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz); - - // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i); - // psPolynomial2DDump (psf->params->data[i]); - } - - psFree (stats); - psFree (x); - psFree (y); - psFree (z); - psFree (dz); - return (true); -} - -// use the model position parameters to construct a realization of the PSF model -// at the object coordinates -pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) { - - // need to define the relationship between the modelFLT and modelPSF ? - - // find function used to set the model parameters - pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type); - - // allocate a new pmModel to hold the PSF version - pmModel *modelPSF = pmModelAlloc (psf->type); - - // set model parameters for this source based on PSF information - modelFromPSFFunc (modelPSF, modelFLT, psf); - - return (modelPSF); -} Index: unk/psphot/src/pmPSF.h =================================================================== --- /trunk/psphot/src/pmPSF.h (revision 5717) +++ (revision ) @@ -1,27 +1,0 @@ - -# ifndef PM_PSF_H -# define PM_PSF_H - -/** pmPSF data structure - * - * - * - */ -typedef struct -{ - pmModelType type; ///< PSF Model in use - psArray *params; ///< Model parameters (psPolynomial2D) - float chisq; ///< PSF goodness statistic - float ApResid; - float dApResid; - float skyBias; - int nPSFstars; ///< number of stars used to measure PSF -} -pmPSF; - -// pmPSF utilities -pmPSF *pmPSFAlloc (pmModelType type); -bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask); -pmModel *pmModelFromPSF (pmModel *model, pmPSF *psf); - -# endif Index: unk/psphot/src/pmPSFtry.c =================================================================== --- /trunk/psphot/src/pmPSFtry.c (revision 5717) +++ (revision ) @@ -1,338 +1,0 @@ -# include -# include "psLibUtils.h" -# include "pmObjects_EAM.h" -# include "psModulesUtils.h" -# include "pmPSF.h" -# include "pmPSFtry.h" -# include "pmModelGroup.h" - -// ******** pmPSFtry functions ************************************************** -// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely -// * fitted version of the model, the pmPSF fit to the fitted model parameters, -// * and the PSF fits to the source. It also includes the statistics from the -// * fits, both the individual sources, and the collection - -// free a pmPSFtry structure -static void pmPSFtryFree (pmPSFtry *test) { - - if (test == NULL) return; - - psFree (test->psf); - psFree (test->sources); - psFree (test->modelFLT); - psFree (test->modelPSF); - psFree (test->metric); - psFree (test->fitMag); - psFree (test->mask); - return; -} - -// allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) { - - pmModelType type; - - pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); - - // XXX probably need to increment ref counter - type = pmModelSetType (modelName); - test->psf = pmPSFAlloc (type); - test->sources = psMemCopy(sources); - test->modelFLT = psArrayAlloc (sources->n); - test->modelPSF = psArrayAlloc (sources->n); - test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); - test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); - test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); - - for (int i = 0; i < test->modelFLT->n; i++) { - test->mask->data.U8[i] = 0; - test->modelFLT->data[i] = NULL; - test->modelPSF->data[i] = NULL; - test->metric->data.F64[i] = 0; - test->fitMag->data.F64[i] = 0; - } - - psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); - return (test); -} - -// build a pmPSFtry for the given model: -// - fit each source with the free-floating model -// - construct the pmPSF from the collection of models -// - fit each source with the PSF-parameter models -// - measure the pmPSF quality metric (dApResid) - -// sources used in for pmPSFtry may be masked by the analysis -// mask values indicate the reason the source was rejected: - -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) -{ - bool status; - float obsMag; - float fitMag; - float x; - float y; - int Nflt = 0; - int Npsf = 0; - - pmPSFtry *try = pmPSFtryAlloc (sources, modelName); - - // stage 1: fit an independent model (freeModel) to all sources - psTimerStart ("fit"); - for (int i = 0; i < try->sources->n; i++) { - - pmSource *source = try->sources->data[i]; - pmModel *model = pmSourceModelGuess (source, try->psf->type); - x = source->peak->x; - y = source->peak->y; - - // set temporary object mask and fit object - // fit model as FLT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - // exclude the poor fits - if (!status) { - try->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; - psFree (model); - continue; - } - try->modelFLT->data[i] = model; - Nflt ++; - } - psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); - psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); - - // make this optional? - // DumpModelFits (try->modelFLT, "modelsFLT.dat"); - - // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (try->psf, try->modelFLT, try->mask); - - // stage 3: refit with fixed shape parameters - psTimerStart ("fit"); - for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - - pmSource *source = try->sources->data[i]; - pmModel *modelFLT = try->modelFLT->data[i]; - - // set shape for this model based on PSF - pmModel *modelPSF = pmModelFromPSF (modelFLT, try->psf); - x = source->peak->x; - y = source->peak->y; - - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, modelPSF, true); - - // skip poor fits - if (!status) { - try->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; - psFree (modelPSF); - goto next_source; - } - - // otherwise, save the resulting model - try->modelPSF->data[i] = modelPSF; - - // XXX : use a different aperture radius from the fit radius? - // XXX : use a different estimator for the local sky? - // XXX : pass 'source' as input? - if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { - try->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; - goto next_source; - } - - try->metric->data.F64[i] = obsMag - fitMag; - try->fitMag->data.F64[i] = fitMag; - Npsf ++; - - next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - } - psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); - psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); - - // make this optional - // DumpModelFits (try->modelPSF, "modelsPSF.dat"); - - // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric_Alt (try, RADIUS); - psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", - modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias); - - return (try); -} - -bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) { - - float dBin; - int nKeep, nSkip; - - // the measured (aperture - fit) magnitudes (dA == try->metric) - // depend on both the true ap-fit (dAo) and the bias in the sky measurement: - // dA = dAo + dsky/flux - // where flux is the flux of the star - // we fit this trend to find the infinite flux aperture correction (dAo), - // the nominal sky bias (dsky), and the error on dAo - // the values of dA are contaminated by stars with close neighbors in the aperture - // we use an outlier rejection to avoid this bias - - // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64); - for (int i = 0; i < try->sources->n; i++) { - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]); - } - - // find min and max of (1/flux): - psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - psVectorStats (stats, rflux, NULL, try->mask, PSFTRY_MASK_ALL); - - // build binned versions of rflux, metric - dBin = (stats->max - stats->min) / 10.0; - psVector *rfBin = psVectorCreate (NULL, stats->min, stats->max, dBin, PS_TYPE_F64); - psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64); - psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8); - psFree (stats); - - psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin); - - // group data in daBin bins, measure lower 50% mean - for (int i = 0; i < daBin->n; i++) { - - psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64); - tmp->n = 0; - - // accumulate data within bin range - for (int j = 0; j < try->sources->n; j++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[j] & PSFTRY_MASK_ALL) continue; - - // skip points with extreme dA values - if (fabs(try->metric->data.F64[j]) > 0.5) continue; - - // skip points outside of this bin - if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue; - if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue; - - tmp->data.F64[tmp->n] = try->metric->data.F64[j]; - tmp->n ++; - } - - // is this a valid point? - maskB->data.U8[i] = 0; - if (tmp->n < 2) { - maskB->data.U8[i] = 1; - psFree (tmp); - continue; - } - - // dA values are contaminated with low outliers - // measure statistics only on upper 50% of points - // this would be easier if we could sort in reverse: - - psVectorSort (tmp, tmp); - nKeep = 0.5*tmp->n; - nSkip = tmp->n - nKeep; - - psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64); - for (int j = 0; j < tmp2->n; j++) { - tmp2->data.F64[j] = tmp->data.F64[j + nSkip]; - } - - stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, tmp2, NULL, NULL, 0); - psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); - - daBin->data.F64[i] = stats->sampleMedian; - - psFree (stats); - psFree (tmp); - psFree (tmp2); - } - - // linear clipped fit to rfBin, daBin - psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); - psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); - - psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); - psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); - - stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); - stats = psVectorStats (stats, daResid, NULL, maskB, 1); - - try->psf->ApResid = poly->coeff[0]; - try->psf->dApResid = stats->clippedStdev; - try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); - - psFree (rflux); - psFree (rfBin); - psFree (daBin); - psFree (maskB); - psFree (daBinFit); - psFree (daResid); - psFree (poly); - psFree (stats); - psFree (fitstat); - - return true; -} - -bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) { - - // the measured (aperture - fit) magnitudes (dA == try->metric) - // depend on both the true ap-fit (dAo) and the bias in the sky measurement: - // dA = dAo + dsky/flux - // where flux is the flux of the star - // we fit this trend to find the infinite flux aperture correction (dAo), - // the nominal sky bias (dsky), and the error on dAo - // the values of dA are contaminated by stars with close neighbors in the aperture - // we use an outlier rejection to avoid this bias - - // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64); - for (int i = 0; i < try->sources->n; i++) { - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]); - } - - // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - // XXX EAM - stats->min = 1.0; - stats->max = 3.0; - stats->clipIter = 3; - - // linear clipped fit to rfBin, daBin - psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); - poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux); - fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - try->psf->ApResid = poly->coeff[0]; - try->psf->dApResid = stats->sampleStdev; - try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); - - FILE *f; - f = fopen ("apresid.dat", "w"); - if (f == NULL) psAbort ("pmPSFtry", "can't open output file"); - - for (int i = 0; i < try->sources->n; i++) { - fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]); - } - fclose (f); - - psFree (rflux); - psFree (poly); - psFree (stats); - - // psFree (daFit); - // psFree (daResid); - - return true; -} Index: unk/psphot/src/pmPSFtry.h =================================================================== --- /trunk/psphot/src/pmPSFtry.h (revision 5717) +++ (revision ) @@ -1,31 +1,0 @@ - -# ifndef PM_PSF_TRY_H -# define PM_PSF_TRY_H - -// data to try a given PSF model type -typedef struct { - 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; -} pmPSFtry; - -enum { - PSFTRY_MASK_CLEAR = 0x00, - PSFTRY_MASK_OUTLIER = 0x01, // 1: outlier in psf polynomial fit (provided by psPolynomials) - PSFTRY_MASK_FLT_FAIL = 0x02, // 2: flt model failed to converge - PSFTRY_MASK_PSF_FAIL = 0x04, // 3: psf model failed to converge - PSFTRY_MASK_BAD_PHOT = 0x08, // 4: invalid source photometry - PSFTRY_MASK_ALL = 0x0f, -} pmPSFtryMaskValues; - -// pmPSFtry utilities -pmPSFtry *pmPSFtryAlloc (psArray *stars, char *modelName); -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float radius); -bool pmPSFtryMetric (pmPSFtry *try, float RADIUS); -bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS); - -# endif Index: /trunk/psphot/src/pmSourceContour.c =================================================================== --- /trunk/psphot/src/pmSourceContour.c (revision 5717) +++ /trunk/psphot/src/pmSourceContour.c (revision 5718) @@ -100,5 +100,5 @@ new implementation of source contour function *****************************************************************************/ -psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold) { +psArray *pmSourceContour_EAM (psImage *image, int xc, int yc, float threshold) { int xR, yR, x0, x1, x0s, x1s; @@ -106,4 +106,7 @@ psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(image, false); + + int x = xc - image->col0; + int y = yc - image->row0; // the requested point must be within the contour @@ -141,5 +144,5 @@ x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); if (x0 == x1) { - fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); goto pt1; } @@ -150,5 +153,5 @@ x1 = findContourNeg (image, threshold, xR, yR, RIGHT); } - fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); xVec->data.F32[Npt + 0] = image->col0 + x0; @@ -175,5 +178,5 @@ x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); if (x0 == x1) { - fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); goto pt2; } @@ -184,5 +187,5 @@ x1 = findContourNeg (image, threshold, xR, yR, RIGHT); } - fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); + // fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); xVec->data.F32[Npt + 0] = image->col0 + x0; @@ -202,5 +205,5 @@ yVec->n = Npt; - fprintf (stderr, "done\n"); + // fprintf (stderr, "done\n"); psArray *tmpArray = psArrayAlloc(2); tmpArray->data[0] = (psPtr *) xVec; Index: unk/psphot/src/psEllipse.c =================================================================== --- /trunk/psphot/src/psEllipse.c (revision 5717) +++ (revision ) @@ -1,59 +1,0 @@ -# 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: unk/psphot/src/psEllipse.h =================================================================== --- /trunk/psphot/src/psEllipse.h (revision 5717) +++ (revision ) @@ -1,24 +1,0 @@ -// 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/psModulesUtils.h =================================================================== --- /trunk/psphot/src/psModulesUtils.h (revision 5717) +++ /trunk/psphot/src/psModulesUtils.h (revision 5718) @@ -5,7 +5,7 @@ # include // for strcasecmp # include +# include +# include # include "psLibUtils.h" -# include "pmObjects_EAM.h" -# include "pmModelGroup.h" // psModule extra utilities Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5717) +++ /trunk/psphot/src/psphot.c (revision 5718) @@ -47,4 +47,7 @@ if (!strcasecmp (breakPt, "CLASS")) exit (0); + psphotBasicDeblend (sources, config, sky); + if (!strcasecmp (breakPt, "DEBLEND")) exit (0); + // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, psphotSortBySN); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5717) +++ /trunk/psphot/src/psphot.h (revision 5718) @@ -21,53 +21,55 @@ // top-level psphot functions -psMetadata *psphotArguments (int *argc, char **argv); -eamReadout *psphotSetup (psMetadata *config); -psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); -psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); -pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); -bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); -void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +psMetadata *psphotArguments (int *argc, char **argv); +eamReadout *psphotSetup (psMetadata *config); +psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); +psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); +void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); -bool psphotSubtractPSF (pmSource *source); -int psphotSortBySN (const void **a, const void **b); -int psphotSortByY (const void **a, const void **b); -int psphotSaveImage (psMetadata *header, psImage *image, char *filename); -bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); +bool psphotSubtractPSF (pmSource *source); +int psphotSortBySN (const void **a, const void **b); +int psphotSortByY (const void **a, const void **b); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +void psphotModelGroupInit (void); -bool pmSourceFitFixed (pmSource *source, pmModel *model); -psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky); -pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); -pmModel *pmSourceSelectModel (pmSource *source); +bool pmSourceFitFixed (pmSource *source, pmModel *model); +psArray *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky); +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +pmModel *pmSourceSelectModel (pmSource *source); // eamReadout functions -eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); +eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); // output functions -bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMP (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteCMF (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); -bool pmPeaksWriteText (psArray *sources, char *filename); -bool pmMomentsWriteText (psArray *sources, char *filename); -bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); -bool pmModelWriteFLTs (psArray *sources, char *filename); -bool pmModelWriteNULLs (psArray *sources, char *filename); +bool pmPeaksWriteText (psArray *sources, char *filename); +bool pmMomentsWriteText (psArray *sources, char *filename); +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); +bool pmModelWriteFLTs (psArray *sources, char *filename); +bool pmModelWriteNULLs (psArray *sources, char *filename); -int pmSourcesDophotType (pmSource *source); +int pmSourcesDophotType (pmSource *source); // psphotModelTest functions -psMetadata *modelTestArguments (int *argc, char **argv); -bool modelTestFitSource (eamReadout *imdata, psMetadata *config); -bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -float psphotCrossProduct (pmSource *Mi, pmSource *Mj); +psMetadata *modelTestArguments (int *argc, char **argv); +bool modelTestFitSource (eamReadout *imdata, psMetadata *config); +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +float psphotCrossProduct (pmSource *Mi, pmSource *Mj); psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); -bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -pmModel *pmModelCopy (pmModel *model); -psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); -psMetadata *psphotTestArguments (int *argc, char **argv); +pmModel *pmModelCopy (pmModel *model); +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); +psMetadata *psphotTestArguments (int *argc, char **argv); +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 5717) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 5718) @@ -1,74 +1,134 @@ # include "psphot.h" -// try PSF models and select best option +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { -pmPSF *psphotBasicDeblend (psMetadata *config, psArray *sources, psStats *skystats) -{ + int N; + bool status; + float threshold; + pmSource *source, *testSource; - // source analysis is done in S/N order (brightest first) + FILE *f = fopen ("deblend.dat", "w"); + if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file"); + int Nblend = 0; + + psTimerStart ("psphot"); + + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + + float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION"); + float NSIGMA = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA"); + float minThreshold = NSIGMA*sky->sampleStdev; + fprintf (stderr, "min threshold: %f\n", minThreshold); + + // we need sources spatially-sorted to find overlaps sources = psArraySort (sources, psphotSortByY); + // source analysis is done in peak order (brightest first) + // we use an index for this so the spatial sorting is kept psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32); for (int i = 0; i < SN->n; i++) { - source = sources->data[i]; - SN->data.F32[i] = source->moments->SN; + source = sources->data[i]; + SN->data.F32[i] = source->peak->counts; } psVector *index = psVectorSortIndex (NULL, SN); - + // this results in an index of increasing SN + + // temporary array for overlapping objects we find psArray *overlap = psArrayAlloc (100); - for (int i = 0; i < sources->n; i++) { - N = index->data.U32[i]; - source = sources->data[N]; - overlap->n = 0; + // XXX make sure this results in decreasing, not increasing, SN + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; - // search backwards for overlapping sources - for (int j = N - 1; j >= 0; j--) { - testSource = sources->data[j]; - if (testSource->moments->x < source->pixels->col0) continue; - if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue; - if (testSource->moments->y < source->pixels->row0) break; - if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) { - fprintf (stderr, "warning: invalid condition\n"); - continue; + if (source->type == PM_SOURCE_BLEND) continue; + + overlap->n = 0; + + // search backwards for overlapping sources + for (int j = N - 1; j >= 0; j--) { + testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) break; + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + psArrayAdd (overlap, 100, testSource); } - psArrayAdd (overlap, 100, testSource); - } - // search forwards for overlapping sources - for (int j = N + 1; j < sources->n; j++) { - testSource = sources->data[j]; - if (testSource->moments->x < source->pixels->col0) continue; - if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue; - if (testSource->moments->y < source->pixels->row0) { - fprintf (stderr, "warning: invalid condition\n"); - continue; + // search forwards for overlapping sources + for (int j = N + 1; j < sources->n; j++) { + testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break; + psArrayAdd (overlap, 100, testSource); } - if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) break; - psArrayAdd (overlap, 100, testSource); - } - // generate source contour (1/4 peak counts) - if (overlap->n > 0) { - // XXX EAM : make the 1/4 user-defined. - threshold = 0.25 * (source->moments->Peak - source->moments->Sky) + source->moments->Sky); - psArray *contour = pmSourceContour_EAM (source->pixels, source->moments->x, source->moments->x, threshold); + // generate source contour (1/4 peak counts) + if (overlap->n > 0) { + // XXX EAM : make the 1/4 user-defined. + // XXX keep threshold from dropping too low (N*sky.sigma) + threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; + threshold = PS_MAX (threshold, minThreshold); - // XXX EAM : the order is wrong: need to check j, j+1 contour points - psVector *xv = contour->data[0]; - psVector *yv = contour->data[1]; - for (int j = 0; j < xv->n; j++) { - for (int k = 0; k < overlap->n; k++) { - testSource = overlap->data[k]; - if (fabs(yv->data.F32[j] - testSource->moments->y) > 0.5) continue; - if (xv->data.F32[j] > testSource->moments->x) continue; - if (xv->data.F32[j] < testSource->moments->x) continue; - // XXX EAM : mark source with flag (blended) - } - } + // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now + psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold); + if (contour == NULL) { + fprintf (stderr, "below threshold? invalid peak?\n"); + continue; + } + + // the source contour consists of two vectors, xv and yv. the contour is + // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]). both + // coordinate pairs have the same yv[] value, with xv[i] corresponding to + // the left boundary and xv[i+1] corresponding to the right boundary + + psVector *xv = contour->data[0]; + psVector *yv = contour->data[1]; + for (int k = 0; k < overlap->n; k++) { + testSource = overlap->data[k]; + if (testSource->peak->counts > source->peak->counts) continue; + for (int j = 0; j < xv->n; j+=2) { + if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; + if (xv->data.F32[j+0] > testSource->peak->x) break; + if (xv->data.F32[j+1] < testSource->peak->x) break; + + int xp0 = source->moments->x - source->pixels->col0; + int xp1 = source->peak->x - source->pixels->col0; + int xp2 = testSource->moments->x - source->pixels->col0; + int xp3 = testSource->peak->x - source->pixels->col0; + + int yp0 = source->moments->y - source->pixels->row0; + int yp1 = source->peak->y - source->pixels->row0; + int yp2 = testSource->moments->y - testSource->pixels->row0; + int yp3 = testSource->peak->y - testSource->pixels->row0; + + fprintf (f, "%d %d (%f, %f) : %d %d (%f, %f) vs %f\n", + source->peak->x, source->peak->y, + source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], + testSource->peak->x, testSource->peak->y, + testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold + ); + + testSource->type = PM_SOURCE_BLEND; + Nblend ++; + j = xv->n; + } + } + } } + fclose (f); + psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot")); + return true; +} - psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); - return (psf); -} + Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 5717) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5718) @@ -13,4 +13,7 @@ sources = psArraySort (sources, psphotSortByY); + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); @@ -23,6 +26,17 @@ // pre-calculate all model pixels psArray *models = psArrayAlloc (sources->n); + psVector *index = psVectorAlloc (sources->n, PS_TYPE_U32); + models->n = 0; + index->n = 0; + for (int i = 0; i < sources->n; i++) { pmSource *inSource = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (inSource->type == PM_SOURCE_BLEND) continue; + if (inSource->type == PM_SOURCE_DEFECT) continue; + if (inSource->type == PM_SOURCE_SATURATED) continue; + pmSource *otSource = pmSourceAlloc (); @@ -72,5 +86,6 @@ // save source in list otSource->modelPSF = model; - models->data[i] = otSource; + index->data.U32[models->n] = i; + psArrayAdd (models, 100, otSource); } psLogMsg ("psphot.emsemble", 4, "built models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n); @@ -79,7 +94,8 @@ // fill out the sparse matrix - psSparse *sparse = psSparseAlloc (sources->n, 100); - for (int i = 0; i < sources->n; i++) { - pmSource *Fi = sources->data[i]; + psSparse *sparse = psSparseAlloc (models->n, 100); + for (int i = 0; i < models->n; i++) { + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; pmSource *Mi = models->data[i]; @@ -93,13 +109,12 @@ // loop over all other stars following this one - for (int j = i + 1; j < sources->n; j++) { + for (int j = i + 1; j < models->n; j++) { pmSource *Mj = models->data[j]; - // XXX double check that this is working! should not break here - check me here; - if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) break; + // skip over disjoint source images, break after last possible overlap + if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; + if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; + if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; - if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) continue; - if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; // got an overlap; calculate cross-product and add to output array @@ -115,6 +130,7 @@ // adjust models, set sources and subtract - for (int i = 0; i < sources->n; i++) { - pmSource *Fi = sources->data[i]; + for (int i = 0; i < models->n; i++) { + int N = index->data.U32[i]; + pmSource *Fi = sources->data[N]; pmSource *Mi = models->data[i]; Index: /trunk/psphot/src/psphotImageBackground.c =================================================================== --- /trunk/psphot/src/psphotImageBackground.c (revision 5717) +++ /trunk/psphot/src/psphotImageBackground.c (revision 5718) @@ -58,7 +58,5 @@ for (int j = 0; j < image->numRows; j++) { for (int i = 0; i < image->numCols; i++) { - if (mask->data.U8[j][i]) { - image->data.F32[j][i] = 0; - } else { + if (!mask->data.U8[j][i]) { Sky = psPolynomial2DEval (skyModel, i, j); image->data.F32[j][i] -= Sky; @@ -66,4 +64,5 @@ } } + psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); Index: /trunk/psphot/src/psphotModelGroupInit.c =================================================================== --- /trunk/psphot/src/psphotModelGroupInit.c (revision 5717) +++ /trunk/psphot/src/psphotModelGroupInit.c (revision 5718) @@ -1,8 +1,8 @@ # include "psphot.h" -# include "models/pmModel_WAUSS.c" +# include "models/pmModel_TAUSS.c" static pmModelGroup userModels[] = { - {"PS_MODEL_WAUSS", 9, pmModelFunc_WAUSS, pmModelFlux_WAUSS, pmModelRadius_WAUSS, pmModelLimits_WAUSS, pmModelGuess_WAUSS, pmModelFromPSF_WAUSS, pmModelFitStatus_WAUSS}, + {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, }; @@ -17,4 +17,4 @@ pmModelGroupAdd (&userModels[i]); } - return true; + return; } Index: /trunk/psphot/src/psphotReapplyPSF.c =================================================================== --- /trunk/psphot/src/psphotReapplyPSF.c (revision 5717) +++ /trunk/psphot/src/psphotReapplyPSF.c (revision 5718) @@ -14,5 +14,8 @@ sources = psArraySort (sources, psphotSortBySN); - // we may set this differently here from the value used to mark likely saturated stars + // this should be added to the PM_SOURCE flags: + int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1; + + // we may set this differently here from the value used to mark likely saturated stars float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); @@ -26,4 +29,5 @@ // skip non-astronomical objects (very likely defects) // XXX EAM : should we try these anyway? + if (source->type == PM_SOURCE_BLEND) continue; if (source->type == PM_SOURCE_DEFECT) continue; if (source->type == PM_SOURCE_SATURATED) continue;