Index: /branches/eam_branch_20070817/psModules/src/camera/pmFPA.h =================================================================== --- /branches/eam_branch_20070817/psModules/src/camera/pmFPA.h (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/camera/pmFPA.h (revision 14546) @@ -6,6 +6,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-14 03:17:18 $ + * @version $Revision: 1.15.4.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -18,4 +18,20 @@ #define FPA_ASTROM 1 ///< Include astrometry information in the structures? + +// Return chip position, given FPA position; calculations are all done in pixel units +#define PM_FPA_TO_CHIP(pos, chip0, chipParity) \ + (((pos) - (chip0))*(chipParity)) + +// Return cell position, given chip position; calculations are all done in pixel units +#define PM_CHIP_TO_CELL(pos, cell0, cellParity, binning) \ + (((pos) - (cell0))*(cellParity)/(binning)) + +// Return chip position, given a cell position; calculations are all done in pixel units +#define PM_CELL_TO_CHIP(pos, cell0, cellParity, binning) \ + ((pos)*(binning)*(cellParity) + (cell0)) + +// Return FPA position, given a chip position; calculations are all done in pixel units +#define PM_CHIP_TO_FPA(pos, chip0, chipParity) \ + ((pos)*(chipParity) + (chip0)) /// Focal plane array (the entirety of the camera) Index: /branches/eam_branch_20070817/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/Makefile.am (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/Makefile.am (revision 14546) @@ -14,4 +14,5 @@ pmSourceContour.c \ pmSourceFitModel.c \ + pmSourceFitSet.c \ pmSourcePhotometry.c \ pmSourceIO.c \ @@ -51,4 +52,5 @@ pmSourceContour.h \ pmSourceFitModel.h \ + pmSourceFitSet.h \ pmSourcePhotometry.h \ pmSourceIO.h \ Index: /branches/eam_branch_20070817/psModules/src/objects/pmModel.h =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmModel.h (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/pmModel.h (revision 14546) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.11.6.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-17 21:01:59 $ + * @version $Revision: 1.11.6.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -169,4 +169,20 @@ ); +bool pmModelAddWithOffset(psImage *image, + psImage *mask, + pmModel *model, + pmModelOpMode mode, + psMaskType maskVal, + int dx, + int dy); + +bool pmModelSubWithOffset(psImage *image, + psImage *mask, + pmModel *model, + pmModelOpMode mode, + psMaskType maskVal, + int dx, + int dy); + /** pmModelFitStatus() * Index: /branches/eam_branch_20070817/psModules/src/objects/pmSource.c =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSource.c (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSource.c (revision 14546) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.34.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-17 21:01:59 $ + * @version $Revision: 1.34.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -797,5 +797,5 @@ // should we call pmSourceCacheModel if it does not exist? -bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal) { +bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy) { bool status; @@ -856,8 +856,9 @@ target = source->weight; } + if (add) { - status = pmModelAdd (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); + status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); } else { - status = pmModelSub (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal); + status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy); } @@ -866,9 +867,17 @@ bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal) { - return pmSourceOp (source, mode, true, maskVal); + return pmSourceOp (source, mode, true, maskVal, 0, 0); } bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal) { - return pmSourceOp (source, mode, false, maskVal); + return pmSourceOp (source, mode, false, maskVal, 0, 0); +} + +bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) { + return pmSourceOp (source, mode, true, maskVal, dx, dy); +} + +bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) { + return pmSourceOp (source, mode, false, maskVal, dx, dy); } Index: /branches/eam_branch_20070817/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSource.h (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSource.h (revision 14546) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.16.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-17 21:01:59 $ + * @version $Revision: 1.16.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -218,6 +218,8 @@ bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal); bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal); - -bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal); +bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy); +bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy); + +bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy); bool pmSourceCacheModel (pmSource *source, psMaskType maskVal); bool pmSourceCachePSF (pmSource *source, psMaskType maskVal); Index: /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitModel.c (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitModel.c (revision 14546) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.24.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-17 21:01:59 $ + * @version $Revision: 1.24.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -230,352 +230,2 @@ } -# define SKIP_FIT_SET 1 -# if (SKIP_FIT_SET) - -bool pmSourceFitSet (pmSource *source, - psArray *modelSet, - pmSourceFitMode mode, - psMaskType maskVal) -{ - psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->maskObj, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - return true; -} - -# else - -/********************* Source Model Set Functions ***************************/ - -// these static variables are used by one pass of pmSourceFitSet to store -// data for a model set. If we are going to make psphot thread-safe, these -// will have to go in a structure of their own or be allocated once per thread -// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, -// nPar = nSrc*(nOnePar - 1) + 1 -static pmModelFunc oneModelFunc; -static pmModelLimits oneCheckLimits; -static psVector *onePar; -static psVector *oneDeriv; -static int nOnePar; - -bool pmSourceFitSetInit (pmModelType type) -{ - - oneModelFunc = pmModelFunc_GetFunction (type); - oneCheckLimits = pmModelLimits_GetFunction (type); - nOnePar = pmModelClassParameterCount (type); - - onePar = psVectorAlloc (nOnePar, PS_DATA_F32); - oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32); - - return true; -} - -void pmSourceFitSetClear (void) -{ - - psFree (onePar); - psFree (oneDeriv); - return; -} - -bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas) -{ - // convert the value of nParam into corresponding single model parameter entry - // convert params into corresponding single model parameter pointer - - int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1; - float *paramSingle = params + nParam - nParamSingle; - float *betaSingle = betas + nParam - nParamSingle; - bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle); - return status; -} - -psF32 pmSourceFitSet_Function(psVector *deriv, - const psVector *params, - const psVector *x) -{ - - psF32 value; - psF32 model; - - psF32 *PAR = onePar->data.F32; - psF32 *dPAR = oneDeriv->data.F32; - - psF32 *pars = params->data.F32; - psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; - - int nSrc = (params->n - 1) / (nOnePar - 1); - - PAR[0] = model = pars[0]; - for (int i = 0; i < nSrc; i++) { - int nOff = i*nOnePar - i; - for (int n = 1; n < nOnePar; n++) { - PAR[n] = pars[nOff + n]; - } - if (deriv == NULL) { - value = oneModelFunc (NULL, onePar, x); - } else { - value = oneModelFunc (oneDeriv, onePar, x); - for (int n = 1; n < nOnePar; n++) { - dpars[nOff + n] = dPAR[n]; - } - } - model += value; - } - if (deriv != NULL) { - dpars[0] = dPAR[0]*2.0; - } - return (model); -} - -/* -i: 0 1 2 -n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 -i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 -*/ - -bool pmSourceFitSet (pmSource *source, - psArray *modelSet, - pmSourceFitMode mode, - psMaskType maskVal) -{ - psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->maskObj, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - - // maximum number of valid pixels - psS32 nPix = source->pixels->numRows * source->pixels->numCols; - - // construct the coordinate and value entries - psArray *x = psArrayAllocEmpty(nPix); - psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32); - psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32); - - // fill in the coordinate and value entries - nPix = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - // skip masked points - if (source->maskObj->data.U8[i][j] & maskVal) { - continue; - } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { - continue; - } - // skip nan values in image - if (!isfinite(source->pixels->data.F32[i][j])) { - continue; - } - - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[nPix] = (psPtr *) coord; - y->data.F32[nPix] = source->pixels->data.F32[i][j]; - - // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky - // as weight to avoid the bias from systematic errors here we would just use the - // source sky variance - if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; - } else { - yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; - } - nPix++; - } - } - x->n = nPix; - y->n = nPix; - yErr->n = nPix; - - // base values on first model (all models must be identical) - pmModel *model = modelSet->data[0]; - - // determine number of model parameters - int nSrc = modelSet->n; - int nPar = model->params->n - 1; // number of object parameters (excluding sky) - - // define parameter vectors for source set - psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - - // set the static variables - pmSourceFitSetInit (model->type); - - // create the minimization constraints - psMinConstraint *constraint = psMinConstraintAlloc(); - constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8); - constraint->checkLimits = pmSourceFitSet_CheckLimits; - - // set the parameter guesses for the multiple models - // first the value for the single sky parameter - params->data.F32[0] = model->params->data.F32[0]; - - // next, the values for the source parameters - for (int i = 0; i < nSrc; i++) { - model = modelSet->data[i]; - for (int n = 1; n < nPar + 1; n++) { - params->data.F32[i*nPar + n] = model->params->data.F32[n]; - } - } - - if (psTraceGetLevel("psModules.objects") >= 5) { - for (int i = 0; i < params->n; i++) { - fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]); - } - } - - // set the parameter masks based on the fitting mode - int nParams = 0; - switch (mode) { - case PM_SOURCE_FIT_NORM: - // NORM-only model fits only source normalization (Io) - nParams = nSrc; - psVectorInit (constraint->paramMask, 1); - for (int i = 0; i < nSrc; i++) { - constraint->paramMask->data.U8[1 + i*nPar] = 0; - } - break; - case PM_SOURCE_FIT_PSF: - // PSF model only fits x,y,Io - nParams = 3*nSrc; - psVectorInit (constraint->paramMask, 1); - for (int i = 0; i < nSrc; i++) { - constraint->paramMask->data.U8[1 + i*nPar] = 0; - constraint->paramMask->data.U8[2 + i*nPar] = 0; - constraint->paramMask->data.U8[3 + i*nPar] = 0; - } - break; - case PM_SOURCE_FIT_EXT: - // EXT model fits all params (except sky) - nParams = nPar*nSrc; - psVectorInit (constraint->paramMask, 0); - constraint->paramMask->data.U8[0] = 1; - break; - default: - psAbort("invalid fitting mode"); - } - - // force the floating parameters to fall within the contraint ranges - for (int i = 0; i < params->n; i++) { - pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); - pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); - } - - if (nPix < nParams + 1) { - psTrace (__func__, 4, "insufficient valid pixels\n"); - psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__); - model->flags |= PM_MODEL_STATUS_BADARGS; - psFree (x); - psFree (y); - psFree (yErr); - psFree (params); - psFree (dparams); - psFree(constraint); - pmSourceFitSetClear (); - return(false); - } - - psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); - - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - - fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function); - if (!fitStatus) { - psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc); - } - - // parameter errors from the covariance matrix - for (int i = 0; i < dparams->n; i++) { - if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); - } - - // get the Gauss-Newton distance for fixed model parameters - if (constraint->paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); - psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8); - altmask->data.U8[0] = 1; - for (int i = 1; i < dparams->n; i++) { - altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1; - } - psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function); - for (int i = 0; i < dparams->n; i++) { - if (!constraint->paramMask->data.U8[i]) - continue; - // note that delta is the value *subtracted* from the parameter - // to get the new guess. for dparams to represent the direction - // of motion, we need to take -delta - dparams->data.F32[i] = -delta->data.F32[i]; - } - psFree (delta); - psFree (altmask); - } - - // assign back the parameters to the models - for (int i = 0; i < nSrc; i++) { - model = modelSet->data[i]; - model->params->data.F32[0] = params->data.F32[0]; - for (int n = 1; n < nPar + 1; n++) { - if (psTraceGetLevel("psModules.objects") >= 4) { - fprintf (stderr, "%f ", params->data.F32[i*nPar + n]); - } - model->params->data.F32[n] = params->data.F32[i*nPar + n]; - model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; - } - psTrace ("psModules.objects", 4, " src %d", i); - - // save the resulting chisq, nDOF, nIter - // these are not unique for any one source - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // set the model success or failure status - model->flags |= PM_MODEL_STATUS_FITTED; - if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; - - // models can go insane: reject these - onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0); - onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->col0 + source->pixels->numCols); - onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0); - onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE; - } - psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); - - source->mode |= PM_SOURCE_MODE_FITTED; - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myMin); - psFree(covar); - psFree(constraint); - psFree(params); - psFree(dparams); - - // free static memory used by pmSourceFitSet - pmSourceFitSetClear (); - - rc = (onPic && fitStatus); - psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); - psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc); - return(rc); -} - -# endif Index: /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.c (revision 14546) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.c (revision 14546) @@ -0,0 +1,417 @@ +/** @file pmSourceFitModel.c + * + * fit single source models to image pixels + * + * @author EAM, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.3.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceFitModel.h" +#include "pmSourceFitSet.h" + +// save as static values so they may be set externally +static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; +static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; +static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; +static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; + +/********************* Source Model Set Functions ***************************/ + +static pmSourceFitSetData *thisSet = NULL; + +static void pmSourceFitSetDataFree (pmSourceFitSetData *set) { + if (!set) return; + + psFree (set->modelSet); + psFree (set->paramSet); + psFree (set->derivSet); + return; +} + +pmSourceFitSetData *pmSourceFitSetDataAlloc (psArray *modelSet) { + + pmSourceFitSetData *set = (pmSourceFitSetData *) psAlloc(sizeof(pmSourceFitSetData)); + psMemSetDeallocator(set, (psFreeFunc) pmSourceFitSetDataFree); + + set->modelSet = psMemIncrRefCounter (modelSet); + set->paramSet = psArrayAlloc (modelSet->n); + set->derivSet = psArrayAlloc (modelSet->n); + set->nParamSet = 0; + + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + + int nParams = pmModelClassParameterCount (model->type); + + set->paramSet->data[i] = psVectorAlloc (nParams, PS_DATA_F32); + set->derivSet->data[i] = psVectorAlloc (nParams, PS_DATA_F32); + + set->nParamSet += nParams; + } + + return set; +} + +// this function is called with the full set of parameters and the beta values in a single vector +bool pmSourceFitSetCheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas) +{ + assert (thisSet); + + // nParam is the parameter in the full sequence. determine which single model this comes from + int nModel = -1; + int nParamOne = -1; + int nParamBase = 0; + for (int i = 0; i < thisSet->modelSet->n; i++) { + psVector *param = thisSet->paramSet->data[i]; + if ((nParamBase <= nParam) && (nParam < nParamBase + param->n)) { + nModel = i; + nParamOne = nParam - nParamBase; + break; + } + nParamBase += param->n; + } + assert (nModel > -1); + + pmModel *model = thisSet->modelSet->data[nModel]; + + // pass the single model function a pointer to the start of that model's sequence + float *paramOne = params + nParamBase; + float *betaOne = betas + nParamBase; + bool status = model->modelLimits (mode, nParamOne, paramOne, betaOne); + return status; +} + +bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set) { + + assert (set); + assert (set->paramSet->n == set->derivSet->n); + + int n = 0; + for (int i = 0; i < set->paramSet->n; i++) { + + psVector *paramOne = set->paramSet->data[i]; + psVector *derivOne = set->derivSet->data[i]; + assert ((deriv == NULL) || (paramOne->n == derivOne->n)); + + for (int j = 0; j < paramOne->n; j++, n++) { + param->data.F32[n] = paramOne->data.F32[j]; + if (deriv) { + deriv->data.F32[n] = derivOne->data.F32[j]; + } + } + } + return true; +} + +bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) { + + assert (set); + assert (set->paramSet->n == set->derivSet->n); + + int n = 0; + for (int i = 0; i < set->paramSet->n; i++) { + + psVector *paramOne = set->paramSet->data[i]; + psVector *derivOne = set->derivSet->data[i]; + assert ((deriv == NULL) || (paramOne->n == derivOne->n)); + + for (int j = 0; j < paramOne->n; j++, n++) { + paramOne->data.F32[j] = param->data.F32[n]; + if (deriv) { + derivOne->data.F32[j] = deriv->data.F32[n]; + } + } + } + return true; +} + +bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus) { + + bool onPic = true; + + assert (set); + + int n = 0; + for (int i = 0; i < set->paramSet->n; i++) { + + pmModel *model = set->modelSet->data[i]; + + for (int j = 0; j < model->params->n; j++, n++) { + if (psTraceGetLevel("psModules.objects") >= 4) { + fprintf (stderr, "%f ", param->data.F32[n]); + } + model->params->data.F32[j] = param->data.F32[n]; + model->dparams->data.F32[j] = dparam->data.F32[n]; + } + psTrace ("psModules.objects", 4, " src %d", i); + + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = nPix - model->params->n; + + // set the model success or failure status + model->flags |= PM_MODEL_STATUS_FITTED; + if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0); + onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0); + onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE; + } + return true; +} + +psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x) +{ + float chisqSum = 0.0; + float chisqOne = 0.0; + + assert (thisSet); + pmSourceFitSetSplit (thisSet, deriv, param); + + for (int i = 0; i < thisSet->modelSet->n; i++) { + + pmModel *model = thisSet->modelSet->data[i]; + + psVector *paramOne = thisSet->paramSet->data[i]; + psVector *derivOne = thisSet->derivSet->data[i]; + + chisqOne = model->modelFunc (derivOne, paramOne, x); + chisqSum += chisqOne; + } + + return (chisqSum); +} + +// XXX allow the mode to be a function of the object (eg, S/N) +bool pmSourceFitSetMasks (psMinConstraint *constraint, pmSourceFitSetData *set, pmSourceFitMode mode) { + + // unmask everyone + psVectorInit (constraint->paramMask, 0); + + int n = 0; + for (int i = 0; i < set->paramSet->n; i++) { + psVector *paramOne = set->paramSet->data[i]; + + switch (mode) { + case PM_SOURCE_FIT_NORM: + // mask all but Xo,Yo,Io + for (int j = 0; j < paramOne->n; j++) { + if (j == PM_PAR_I0) continue; + constraint->paramMask->data.U8[n + j] = 1; + } + break; + case PM_SOURCE_FIT_PSF: + // mask all but Xo,Yo,Io + for (int j = 0; j < paramOne->n; j++) { + if (j == PM_PAR_XPOS) continue; + if (j == PM_PAR_YPOS) continue; + if (j == PM_PAR_I0) continue; + constraint->paramMask->data.U8[n + j] = 1; + } + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + constraint->paramMask->data.U8[n + PM_PAR_SKY] = 1; + break; + default: + psAbort("invalid fitting mode"); + } + n += paramOne->n; + } + return true; +} + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + pmSourceFitMode mode, + psMaskType maskVal) +{ + psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->maskObj, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + bool fitStatus = true; + bool onPic = true; + + // maximum number of valid pixels + int nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAllocEmpty(nPix); + psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32); + + // fill in the coordinate and value entries + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->maskObj->data.U8[i][j] & maskVal) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + // skip nan values in image + if (!isfinite(source->pixels->data.F32[i][j])) { + continue; + } + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + + // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky + // as weight to avoid the bias from systematic errors here we would just use the + // source sky variance + if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; + } + nPix++; + } + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + thisSet = pmSourceFitSetDataAlloc (modelSet); + + // define param and deriv vectors for complete set of parameters + psVector *params = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32); + + // set the param and deriv vectors based on the curent values + pmSourceFitSetJoin (NULL, params, thisSet); + + // create the minimization constraints + psMinConstraint *constraint = psMinConstraintAlloc(); + constraint->paramMask = psVectorAlloc (thisSet->nParamSet, PS_TYPE_U8); + constraint->checkLimits = pmSourceFitSetCheckLimits; + + pmSourceFitSetMasks (constraint, thisSet, mode); + + // force the floating parameters to fall within the contraint ranges + for (int i = 0; i < params->n; i++) { + pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); + pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); + } + + if (psTraceGetLevel("psModules.objects") >= 5) { + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]); + } + } + + if (nPix < thisSet->nParamSet + 1) { + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__); + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + model->flags |= PM_MODEL_STATUS_BADARGS; + } + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree(constraint); + psFree (thisSet); + thisSet = NULL; + return(false); + } + + psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); + + fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSetFunction); + if (!fitStatus) { + psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n); + } + + // parameter errors from the covariance matrix + psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32); + for (int i = 0; i < dparams->n; i++) { + if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (constraint->paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); + psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8); + altmask->data.U8[0] = 1; + for (int i = 1; i < dparams->n; i++) { + altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1; + } + psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction); + for (int i = 0; i < dparams->n; i++) { + if (!constraint->paramMask->data.U8[i]) + continue; + // note that delta is the value *subtracted* from the parameter + // to get the new guess. for dparams to represent the direction + // of motion, we need to take -delta + dparams->data.F32[i] = -delta->data.F32[i]; + } + psFree (delta); + psFree (altmask); + } + + pmSourceFitSetValues (thisSet, dparams, params, source, myMin, y->n, fitStatus); + psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n); + + source->mode |= PM_SOURCE_MODE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constraint); + psFree(params); + psFree(dparams); + psFree(thisSet); + + thisSet = NULL; + + bool rc = (onPic && fitStatus); + psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc); + return(rc); +} Index: /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.h =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.h (revision 14546) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSourceFitSet.h (revision 14546) @@ -0,0 +1,53 @@ +/* @file pmSourceFitSet.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.3.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +# ifndef PM_SOURCE_FIT_SET_H +# define PM_SOURCE_FIT_SET_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +typedef struct { + psArray *modelSet; + psArray *paramSet; + psArray *derivSet; + int nParamSet; +} pmSourceFitSetData; + +// initialize data for a group of object models +pmSourceFitSetData *pmSourceFitSetDataAlloc (psArray *modelSet); + +// function used to set limits for a group of models +bool pmSourceFitSetCheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas); + +bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set); +bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param); +bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus); + +psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x); +bool pmSourceFitSetMasks (psMinConstraint *constraint, pmSourceFitSetData *set, pmSourceFitMode mode); + +/** pmSourceFitSet() + * + * Fit the requested model to the specified source. The starting guess for the model is given + * by the input source.model parameter values. The pixels of interest are specified by the + * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image + * data. The function returns TRUE on success or FALSE on failure. + * + */ +bool pmSourceFitSet( + pmSource *source, ///< The input pmSource + psArray *modelSet, ///< model to be fitted + pmSourceFitMode mode, ///< define parameters to be fitted + psMaskType maskVal ///< Vale to mask + +); + +/// @} +# endif /* PM_SOURCE_FIT_MODEL_H */ Index: /branches/eam_branch_20070817/psModules/src/objects/pmSourceUtils.c =================================================================== --- /branches/eam_branch_20070817/psModules/src/objects/pmSourceUtils.c (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/objects/pmSourceUtils.c (revision 14546) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-17 21:01:59 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-20 01:58:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -76,16 +76,32 @@ } - // XXX I don't yet know how this function should be defined -# if (0) - source->peak = pmPeakAlloc (); + pmCell *cell = readout->parent; - float x = model->params->data.F32[PM_PAR_XPOS]; - float y = model->params->data.F32[PM_PAR_YPOS]; + float Io = model->params->data.F32[PM_PAR_I0]; + float xChip = model->params->data.F32[PM_PAR_XPOS]; + float yChip = model->params->data.F32[PM_PAR_YPOS]; - // XXX need to define the radius in some rational way - // XXX x,y are defined wrt readout->image parent, but the model - // parameters are defined wrt chip coordinates - pmSourceDefinePixels (source, readout, x, y, radius); -# endif + source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE); + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + // XXX whose binning? + int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + + // Position on the cell + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); + + // Position on the readout + // float xReadout = CELL_TO_READOUT(xCell, x0Readout); + // float yReadout = CELL_TO_READOUT(yCell, y0Readout); + + float radius = 5.0; // XXX need to define the radius in some rational way + + pmSourceDefinePixels (source, readout, xCell, yCell, radius); return (source); Index: /branches/eam_branch_20070817/psModules/src/psmodules.h =================================================================== --- /branches/eam_branch_20070817/psModules/src/psmodules.h (revision 14545) +++ /branches/eam_branch_20070817/psModules/src/psmodules.h (revision 14546) @@ -88,4 +88,5 @@ #include #include +#include #include #include