Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 14337) +++ trunk/psphot/src/psphot.h (revision 14338) @@ -16,5 +16,5 @@ psString psphotVersionLong(void); -bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark); +bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark); bool psphotReadout (pmConfig *config, const pmFPAview *view); bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources); @@ -116,5 +116,5 @@ bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal); bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal); -psKernel *psphotKernelFromPSF (pmSource *source); +psKernel *psphotKernelFromPSF (pmSource *source, int nPix); bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal); @@ -122,4 +122,15 @@ bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal); bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal); + +// structures & functions to support psf-convolved model fitting + +// pmPCMData : PSF Convolved Model data storage structure +typedef struct { + psImage *model; + psArray *dmodels; + psImage *modelConv; + psArray *dmodelsConv; +} pmPCMData; + // psf-convolved model fitting @@ -129,5 +140,5 @@ psVector *params, psMinConstraint *constraint, - const pmSource *source, + pmSource *source, const psKernel *psf, psMinimizeLMChi2Func func); @@ -138,8 +149,15 @@ const psVector *params, const psVector *paramMask, + pmPCMData *pcm, const pmSource *source, const psKernel *psf, psMinimizeLMChi2Func func); +pmPCMData *pmPCMDataAlloc ( + const psVector *params, + const psVector *paramMask, + pmSource *source); + +psImage *pmPCMDataSaveImage (pmPCMData *pcm); #endif Index: trunk/psphot/src/psphotKernelFromPSF.c =================================================================== --- trunk/psphot/src/psphotKernelFromPSF.c (revision 14337) +++ trunk/psphot/src/psphotKernelFromPSF.c (revision 14338) @@ -1,5 +1,5 @@ # include "psphot.h" -psKernel *psphotKernelFromPSF (pmSource *source) { +psKernel *psphotKernelFromPSF (pmSource *source, int nPix) { assert (source); @@ -10,5 +10,5 @@ // need to decide on the size: dynamically? statically? - psKernel *psf = psKernelAlloc (-1, +1, -1, +1); + psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); for (int j = psf->yMin; j <= psf->yMax; j++) { Index: trunk/psphot/src/psphotModelTest.c =================================================================== --- trunk/psphot/src/psphotModelTest.c (revision 14337) +++ trunk/psphot/src/psphotModelTest.c (revision 14338) @@ -1,11 +1,10 @@ # include "psphotInternal.h" -static char DEFAULT_MODE[] = "EXT"; - -// XXX consider this function : add more test information? -bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { +# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY + +// XXX add more test information? +bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { bool status; int modelType; - unsigned int Nfail; float obsMag, fitMag, value; char name[64]; @@ -13,5 +12,10 @@ pmSourceFitMode fitMode; - psMetadataItem *item = NULL; + // run model fitting tests on a single source? + if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false; + + // find the currently selected readout + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); + PS_ASSERT_PTR_NON_NULL (readout, false); // use poissonian errors or local-sky errors @@ -20,24 +24,41 @@ pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS); - // run model fitting tests on a single source - if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false; + // find the various fitting parameters (try test values first) + float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS"); + if (!status || !isfinite(INNER)) { + INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + } + float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS"); + if (!status || !isfinite(OUTER)) { + OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + } + float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS"); + if (!status || !isfinite(RADIUS)) { + RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + } + float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS"); + if (!status || !isfinite(mRADIUS)) { + mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + } + + // define the source of interest + float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); + float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); + if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined"); // what fitting mode to use? + fitMode = PM_SOURCE_FIT_EXT; char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE"); - if (!status || !strcasecmp (fitModeWord, "DEFAULT")) { - fitModeWord = DEFAULT_MODE; - } - fitMode = PM_SOURCE_FIT_EXT; - if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF; + if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF; + if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT; + if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT; // in fitMode, psf sets the model type if (fitMode == PM_SOURCE_FIT_PSF) { - // XXX load psf using psphotLoadPSF - char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE"); - if (!status) psAbort("PSF_INPUT_FILE not supplied"); - psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE); - psf = pmPSFfromMetadata (psfData); + psf = psphotLoadPSF (config, view, recipe); + if (!psf) psAbort("PSF_INPUT_FILE not supplied"); modelType = psf->type; - } else { + } + if (fitMode == PM_SOURCE_FIT_EXT) { // find the model: supplied by user or first in the PSF_MODEL list char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); @@ -57,5 +78,5 @@ // take the first list element - item = psListGet (list, PS_LIST_HEAD); + psMetadataItem *item = psListGet (list, PS_LIST_HEAD); modelName = item->data.V; } @@ -63,29 +84,32 @@ if (modelType < 0) psAbort("unknown model %s", modelName); } - - // find the fitting parameters (try test values first) - float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS"); - if (!status || !isfinite(INNER)) { - INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - } - - float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS"); - if (!status || !isfinite(OUTER)) { - OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); - } - - float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS"); - if (!status || !isfinite(RADIUS)) { - RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); - } - - float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS"); - if (!status || !isfinite(mRADIUS)) { - mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); - } - - // define the source of interest - float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); - float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); + if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) { + // we need to load BOTH a psf and an ext model + psf = psphotLoadPSF (config, view, recipe); + if (!psf) psAbort("PSF_INPUT_FILE not supplied"); + + // find the model: supplied by user or first in the PSF_MODEL list + char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); + if (!status || !strcasecmp (modelName, "DEFAULT")) { + // get the list pointers for the PSF_MODEL entries + + psList *list = NULL; + psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL"); + if (mdi == NULL) psAbort("missing PSF_MODEL selection"); + if (mdi->type == PS_DATA_STRING) { + list = psListAlloc(NULL); + psListAdd (list, PS_LIST_HEAD, mdi); + } else { + if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection"); + list = psMemIncrRefCounter(mdi->data.list); + } + + // take the first list element + psMetadataItem *item = psListGet (list, PS_LIST_HEAD); + modelName = item->data.V; + } + modelType = pmModelSetType (modelName); + if (modelType < 0) psAbort("unknown model %s", modelName); + } // construct the source structures @@ -100,5 +124,5 @@ // get the source moments status = pmSourceMoments (source, mRADIUS); - if (!status) psAbort("pmSourceLocalSky error"); + if (!status) psAbort("psSourceMoments error"); source->peak->value = source->moments->Peak; @@ -116,16 +140,14 @@ // get the initial model parameter guess pmModel *model = pmSourceModelGuess (source, modelType); - // if any parameters are defined, use those values + + // if any parameters are defined by the user, take those values int nParams = pmModelParameterCount (modelType); psF32 *params = model->params->data.F32; + params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value, + params[PM_PAR_YPOS] = yObj; // XXX or use the centroid for (int i = 0; i < nParams; i++) { - if (i == 2) { - params[i] = xObj; - continue; - } - if (i == 3) { - params[i] = yObj; - continue; - } + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + sprintf (name, "TEST_FIT_PAR%d", i); value = psMetadataLookupF32 (&status, recipe, name); @@ -138,4 +160,5 @@ fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y); + // for PSF fitting, set the shape parameters based on the PSF & source position if (fitMode == PM_SOURCE_FIT_PSF) { pmModel *modelPSF = pmModelFromPSF (model, psf); @@ -154,5 +177,4 @@ fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor); - fprintf (stderr, "input parameters: \n"); for (int i = 0; i < nParams; i++) { @@ -162,4 +184,5 @@ // define the pixels used for the fit psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", mark); + psphotSaveImage (NULL, source->maskObj, "mask1.fits"); char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); @@ -169,5 +192,14 @@ } - status = pmSourceFitModel (source, model, fitMode, maskVal); + if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) { + // build the psf for the object + source->modelPSF = pmModelFromPSF (model, psf); + source->modelEXT = model; + status = psphotPSFConvModel (source, recipe, maskVal); + model = source->modelConv; + params = model->params->data.F32; + } else { + status = pmSourceFitModel (source, model, fitMode, maskVal); + } // measure the source mags Index: trunk/psphot/src/psphotModelWithPSF.c =================================================================== --- trunk/psphot/src/psphotModelWithPSF.c (revision 14337) +++ trunk/psphot/src/psphotModelWithPSF.c (revision 14338) @@ -6,9 +6,9 @@ psVector *params, psMinConstraint *constraint, - const pmSource *source, + pmSource *source, const psKernel *psf, psMinimizeLMChi2Func func) { - psTrace("psLib.math", 3, "---- begin ----\n"); + psTrace("psphot", 3, "---- begin ----\n"); PS_ASSERT_PTR_NON_NULL(min, false); PS_ASSERT_VECTOR_NON_NULL(params, false); @@ -46,6 +46,9 @@ psF32 dLinear = 0.0; + // generate PCM data storage structure + pmPCMData *pcm = pmPCMDataAlloc (params, paramMask, source); + // calculate initial alpha and beta, set chisq (min->value) - min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, source, psf, func); + min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, pcm, source, psf, func); if (isnan(min->value)) { min->iter = min->maxIter; @@ -53,9 +56,9 @@ } // dump some useful info if trace is defined - if (psTraceGetLevel("psLib.math") >= 6) { + if (psTraceGetLevel("psphot") >= 6) { p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)"); p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)"); } - if (psTraceGetLevel("psLib.math") >= 5) { + if (psTraceGetLevel("psphot") >= 5) { p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)"); } @@ -63,6 +66,6 @@ // iterate until the tolerance is reached, or give up while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) { - psTrace("psLib.math", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); - psTrace("psLib.math", 5, "Last delta is %f. Min->tol is %f.\n", min->lastDelta, min->tol); + psTrace("psphot", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); + psTrace("psphot", 5, "Last delta is %f. Min->tol is %f.\n", min->lastDelta, min->tol); @@ -75,14 +78,14 @@ // dump some useful info if trace is defined - if (psTraceGetLevel("psLib.math") >= 6) { + if (psTraceGetLevel("psphot") >= 6) { p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)"); p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)"); } - if (psTraceGetLevel("psLib.math") >= 5) { + if (psTraceGetLevel("psphot") >= 5) { p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)"); } // calculate Chisq for new guess, update Alpha & Beta - Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, source, psf, func); + Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, pcm, source, psf, func); if (isnan(Chisq)) { min->iter ++; @@ -97,9 +100,9 @@ psF32 rho = (min->value - Chisq) / dLinear; - psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, + psTrace("psphot", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho); // dump some useful info if trace is defined - if (psTraceGetLevel("psLib.math") >= 6) { + if (psTraceGetLevel("psphot") >= 6) { p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)"); p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)"); @@ -114,4 +117,8 @@ params = psVectorCopy(params, Params, PS_TYPE_F32); lambda *= 0.1; + + // save the new convolved model image + psFree (source->modelFlux); + source->modelFlux = pmPCMDataSaveImage(pcm); } else { lambda *= 10.0; @@ -119,10 +126,10 @@ min->iter++; } - psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); + psTrace("psphot", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); // construct & return the covariance matrix (if requested) if (covar != NULL) { if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { - psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n"); + psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); } } @@ -134,11 +141,12 @@ psFree(Beta); psFree(Params); + psFree(pcm); if (min->iter == min->maxIter) { - psTrace("psLib.math", 3, "---- end (false) ----\n"); + psTrace("psphot", 3, "---- end (false) ----\n"); return(false); } - psTrace("psLib.math", 3, "---- end (true) ----\n"); + psTrace("psphot", 3, "---- end (true) ----\n"); return(true); } @@ -149,4 +157,5 @@ const psVector *params, const psVector *paramMask, + pmPCMData *pcm, const pmSource *source, const psKernel *psf, @@ -168,14 +177,8 @@ } + // generate the model and derivative images for this parameter set + + // storage for model derivatives psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32); - - // generate the model and derivative images for this parameter set - psImage *model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); - psArray *dmodels = psArrayAlloc (params->n); - for (psS32 n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - psImage *dmodel = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); - dmodels->data[n] = dmodel; - } // working vector to store local coordinate @@ -206,24 +209,34 @@ coord->data.F32[1] = (psF32) (i + source->pixels->row0); - model->data.F32[i][j] = func (deriv, params, coord); + pcm->model->data.F32[i][j] = func (deriv, params, coord); for (int n = 0; n < params->n; n++) { if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - psImage *dmodel = dmodels->data[n]; + psImage *dmodel = pcm->dmodels->data[n]; dmodel->data.F32[i][j] = deriv->data.F32[n]; } } } - psFree (coord); + psFree(coord); + psFree(deriv); + + psphotSaveImage (NULL, pcm->model, "model1.fits"); // convolve model and dmodel arrays with PSF - psImage *modelConv = psImageConvolveDirect (model, psf); - psArray *dmodelsConv = psArrayAlloc (params->n); - for (int n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - psImage *dmodel = dmodels->data[n]; - psImage *dmodelConv = psImageConvolveDirect (dmodel, psf); - dmodelsConv->data[n] = dmodelConv; - } + psImageConvolveDirect (pcm->modelConv, pcm->model, psf); + for (int n = 0; n < pcm->dmodels->n; n++) { + if (pcm->dmodels->data[n] == NULL) continue; + psImage *dmodel = pcm->dmodels->data[n]; + psImage *dmodelConv = pcm->dmodelsConv->data[n]; + psImageConvolveDirect (dmodelConv, dmodel, psf); + } + + // XXX TEST : SAVE IMAGES + psphotSaveImage (NULL, pcm->model, "model.fits"); + psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits"); + psphotSaveImage (NULL, source->pixels, "obj.fits"); + psphotSaveImage (NULL, source->maskObj, "mask.fits"); + psphotSaveImage (NULL, source->weight, "weight.fits"); + exit (0); // zero alpha and beta for summing below @@ -248,5 +261,5 @@ } - float ymodel = modelConv->data.F32[i][j]; + float ymodel = pcm->modelConv->data.F32[i][j]; float yweight = 1.0 / source->weight->data.F32[i][j]; float delta = ymodel - source->pixels->data.F32[i][j]; @@ -263,5 +276,5 @@ continue; } - psImage *dmodel = dmodelsConv->data[n1]; + psImage *dmodel = pcm->dmodelsConv->data[n1]; float weight = dmodel->data.F32[i][j] * yweight; for (psS32 n2 = 0; n2 <= n1; n2++) { @@ -269,5 +282,5 @@ continue; } - dmodel = dmodelsConv->data[n2]; + dmodel = pcm->dmodelsConv->data[n2]; alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j]; } @@ -294,14 +307,50 @@ } - psFree (model); - psFree (dmodels); - psFree (modelConv); - psFree (dmodelsConv); - psFree(deriv); - return(chisq); } - +static void pmPCMDataFree (pmPCMData *pcm) { + + if (pcm == NULL) return; + + psFree (pcm->model); + psFree (pcm->modelConv); + psFree (pcm->dmodels); + psFree (pcm->dmodelsConv); + return; +} + +pmPCMData *pmPCMDataAlloc ( + const psVector *params, + const psVector *paramMask, + pmSource *source) { + + pmPCMData *pcm = (pmPCMData *) psAlloc(sizeof(pmPCMData)); + psMemSetDeallocator(pcm, (psFreeFunc) pmPCMDataFree); + + // Allocate storage images for raw model and derivative images + pcm->model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + pcm->dmodels = psArrayAlloc (params->n); + for (psS32 n = 0; n < params->n; n++) { + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodels->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + } + + // Allocate storage images for convolved model and derivative images + pcm->modelConv = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + pcm->dmodelsConv = psArrayAlloc (params->n); + for (psS32 n = 0; n < params->n; n++) { + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodelsConv->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + } + + return pcm; +} + +psImage *pmPCMDataSaveImage (pmPCMData *pcm) { + + psImage *model = psImageCopy (NULL, pcm->model, PS_TYPE_F32); + return model; +} /* Index: trunk/psphot/src/psphotPSFConvModel.c =================================================================== --- trunk/psphot/src/psphotPSFConvModel.c (revision 14337) +++ trunk/psphot/src/psphotPSFConvModel.c (revision 14338) @@ -7,10 +7,15 @@ // static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; +// input source has both modelPSF and modelEXT. on successful exit, we set the +// modelConv to contain the fitted parameters, and the modelFlux to contain the +// convolved model image. bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) { - // XXX make sure we save a cached copy of the source psf + // make sure we save a cached copy of the psf flux + pmSourceCachePSF (source, maskVal); // convert the cached cached psf model for this source to a psKernel - psKernel *psf = psphotKernelFromPSF (source); + // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius) + psKernel *psf = psphotKernelFromPSF (source, 2); // generate copy of the model @@ -48,4 +53,5 @@ } + // set up the minimization process psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); @@ -54,5 +60,5 @@ bool fitStatus = psphotModelWithPSF_LMM (myMin, covar, params, constraint, source, psf, modelFunc); for (int i = 0; i < dparams->n; i++) { - if (psTraceGetLevel("psModules.objects") >= 4) { + if (psTraceGetLevel("psphot") >= 4) { fprintf (stderr, "%f ", params->data.F32[i]); } @@ -61,5 +67,5 @@ dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); } - psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); + psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); // save the resulting chisq, nDOF, nIter @@ -91,5 +97,5 @@ bool retval = (onPic && fitStatus); - psTrace("psModules.objects", 5, "---- %s(%d) end ----\n", __func__, retval); + psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval); return(retval); }