Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am (revision 28692) @@ -61,4 +61,7 @@ pmPSFtryFitPSF.c \ pmPSFtryMetric.c \ + pmPCMdata.c \ + pmPCM_MinimizeChisq.c \ + pmSourceFitPCM.c \ pmTrend2D.c \ pmGrowthCurveGenerate.c \ @@ -105,4 +108,5 @@ pmPSF_IO.h \ pmPSFtry.h \ + pmPCMdata.h \ pmTrend2D.h \ pmGrowthCurve.h \ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 28692) @@ -1,34 +1,66 @@ -# include "psphotInternal.h" -# define SAVE_IMAGES 0 - -bool psphotModelWithPSF_LMM ( +/* @file pmPCM_MinimizeChisq.c + * structures and functions to support PSF-convolved model fitting + * + * @author EAM, IfA + * + * @version $Revision: 1.29 $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2010 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPCMdata.h" + +# define FACILITY "psModules.objects" + +bool pmPCM_MinimizeChisq ( psMinimization *min, psImage *covar, psVector *params, - psMinConstraint *constraint, pmSource *source, - const psKernel *psf, - psMinimizeLMChi2Func func) + pmPCMdata *pcm) { - psTrace("psphot", 3, "---- begin ----\n"); + psTrace(FACILITY, 3, "---- begin ----\n"); PS_ASSERT_PTR_NON_NULL(min, false); PS_ASSERT_VECTOR_NON_NULL(params, false); PS_ASSERT_VECTOR_NON_EMPTY(params, false); PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); - psVector *paramMask = NULL; - if (constraint != NULL) { - paramMask = constraint->paramMask; - if (paramMask != NULL) { - PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false); - } - } - PS_ASSERT_PTR_NON_NULL(func, false); PS_ASSERT_PTR_NON_NULL(source, false); - - psMinimizeLMLimitFunc checkLimits = NULL; - if (constraint) { - checkLimits = constraint->checkLimits; - } + PS_ASSERT_PTR_NON_NULL(pcm, false); + PS_ASSERT_VECTOR_TYPE(pcm->constraint->paramMask, PS_TYPE_VECTOR_MASK, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(params, pcm->constraint->paramMask, false); + + psVector *paramMask = pcm->constraint->paramMask; + + psMinimizeLMLimitFunc checkLimits = pcm->constraint->checkLimits; // this function has test values and current values for several things @@ -54,9 +86,6 @@ 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, pcm, source, psf, func); + min->value = pmPCM_SetABX(alpha, beta, params, paramMask, pcm, source); if (isnan(min->value)) { min->iter = min->maxIter; @@ -64,21 +93,22 @@ } // dump some useful info if trace is defined - if (psTraceGetLevel("psphot") >= 6) { + if (psTraceGetLevel(FACILITY) >= 6) { p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)"); p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)"); } - if (psTraceGetLevel("psphot") >= 5) { + if (psTraceGetLevel(FACILITY) >= 5) { p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)"); } // iterate until the tolerance is reached, or give up - while ((min->iter < min->maxIter) && ((min->lastDelta > min->minTol) || !isfinite(min->lastDelta))) { - psTrace("psphot", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); - psTrace("psphot", 5, "Last delta is %f. Min->minTol is %f.\n", min->lastDelta, min->minTol); - + bool done = (min->iter >= min->maxIter); + while (!done) { + psTrace("psModules.objects", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); + psTrace("psModules.objects", 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol); // set a new guess for Alpha, Beta, Params if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) { min->iter ++; + if (min->iter >= min->maxIter) break; lambda *= 10.0; continue; @@ -86,17 +116,18 @@ // dump some useful info if trace is defined - if (psTraceGetLevel("psphot") >= 6) { + if (psTraceGetLevel(FACILITY) >= 6) { p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)"); p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)"); p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)"); } - if (psTraceGetLevel("psphot") >= 5) { + if (psTraceGetLevel(FACILITY) >= 5) { p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)"); } // calculate Chisq for new guess, update Alpha & Beta - Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, pcm, source, psf, func); + Chisq = pmPCM_SetABX(Alpha, Beta, Params, paramMask, pcm, source); if (isnan(Chisq)) { min->iter ++; + if (min->iter >= min->maxIter) break; lambda *= 10.0; continue; @@ -109,9 +140,8 @@ psF32 rho = (min->value - Chisq) / dLinear; - psTrace("psphot", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, - Chisq, min->lastDelta, rho); + psTrace(FACILITY, 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("psphot") >= 6) { + if (psTraceGetLevel(FACILITY) >= 6) { p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)"); p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)"); @@ -119,5 +149,5 @@ /* if (Chisq < min->value) { */ - if (rho > 0.0) { + if (rho >= -1e-6) { min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n); min->value = Chisq; @@ -129,16 +159,24 @@ // save the new convolved model image psFree (source->modelFlux); - source->modelFlux = pmPCMDataSaveImage(pcm); + source->modelFlux = pmPCMdataSaveImage(pcm); } else { lambda *= 10.0; } min->iter++; - } - psTrace("psphot", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); + + done = (min->iter >= min->maxIter); + + // check for convergence: + float chisqDOF = Chisq / pcm->nDOF; + if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) { + done |= (min->lastDelta < min->minTol); + } + } + psTrace(FACILITY, 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(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { - psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); + psTrace (FACILITY, 5, "failure to calculate covariance matrix\n"); } // set covar values which are not masked @@ -164,24 +202,21 @@ psFree(Beta); psFree(Params); - psFree(pcm); // if the last improvement was at least as good as maxTol, accept the fit: if (min->lastDelta <= min->maxTol) { - psTrace("psphot", 6, "---- end (true) ----\n"); + psTrace(FACILITY, 6, "---- end (true) ----\n"); return(true); } - psTrace("psphot", 6, "---- end (false) ----\n"); + psTrace(FACILITY, 6, "---- end (false) ----\n"); return(false); } -psF32 psphotModelWithPSF_SetABX( +psF32 pmPCM_SetABX( psImage *alpha, psVector *beta, const psVector *params, const psVector *paramMask, - pmPCMData *pcm, - const pmSource *source, - const psKernel *psf, - psMinimizeLMChi2Func func) + pmPCMdata *pcm, + const pmSource *source) { // XXX: Check vector sizes. @@ -208,8 +243,8 @@ psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - psImageInit (pcm->model, 0.0); + psImageInit (pcm->modelFlux, 0.0); for (int n = 0; n < params->n; n++) { - if (!pcm->dmodels->data[n]) continue; - psImageInit (pcm->dmodels->data[n], 0.0); + if (!pcm->dmodelsFlux->data[n]) continue; + psImageInit (pcm->dmodelsFlux->data[n], 0.0); } @@ -243,9 +278,9 @@ coord->data.F32[1] = (psF32) (i + source->pixels->row0); - pcm->model->data.F32[i][j] = func (deriv, params, coord); + pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord); for (int n = 0; n < params->n; n++) { if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - psImage *dmodel = pcm->dmodels->data[n]; + psImage *dmodel = pcm->dmodelsFlux->data[n]; dmodel->data.F32[i][j] = deriv->data.F32[n]; } @@ -256,17 +291,17 @@ // convolve model and dmodel arrays with PSF - 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); + psImageConvolveDirect (pcm->modelConvFlux, pcm->modelFlux, pcm->psf); + for (int n = 0; n < pcm->dmodelsFlux->n; n++) { + if (pcm->dmodelsFlux->data[n] == NULL) continue; + psImage *dmodel = pcm->dmodelsFlux->data[n]; + psImage *dmodelConv = pcm->dmodelsConvFlux->data[n]; + psImageConvolveDirect (dmodelConv, dmodel, pcm->psf); } // XXX TEST : SAVE IMAGES # if (SAVE_IMAGES) - psphotSaveImage (NULL, psf->image, "psf.fits"); - psphotSaveImage (NULL, pcm->model, "model.fits"); - psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits"); + psphotSaveImage (NULL, pcm->psf->image, "psf.fits"); + psphotSaveImage (NULL, pcm->modelFlux, "model.fits"); + psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits"); psphotSaveImage (NULL, source->pixels, "obj.fits"); psphotSaveImage (NULL, source->maskObj, "mask.fits"); @@ -297,5 +332,5 @@ } - float ymodel = pcm->modelConv->data.F32[i][j]; + float ymodel = pcm->modelConvFlux->data.F32[i][j]; float yweight = 1.0 / source->variance->data.F32[i][j]; float delta = ymodel - source->pixels->data.F32[i][j]; @@ -309,9 +344,9 @@ for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue; - psImage *dmodel = pcm->dmodelsConv->data[n1]; + psImage *dmodel = pcm->dmodelsConvFlux->data[n1]; float weight = dmodel->data.F32[i][j] * yweight; for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue; - dmodel = pcm->dmodelsConv->data[n2]; + dmodel = pcm->dmodelsConvFlux->data[n2]; alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; N2++; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c (revision 28692) @@ -1,47 +1,258 @@ -static void pmPCMDataFree (pmPCMData *pcm) { +/* @file pmPCMdata.c + * structures and functions to support PSF-convolved model fitting + * + * @author EAM, IfA + * + * @version $Revision: 1.29 $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2010 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPCMdata.h" + +static void pmPCMdataFree (pmPCMdata *pcm) { if (pcm == NULL) return; - psFree (pcm->model); + psFree (pcm->modelFlux); + psFree (pcm->modelConvFlux); + psFree (pcm->dmodelsFlux); + psFree (pcm->dmodelsConvFlux); + psFree (pcm->modelConv); - psFree (pcm->dmodels); - psFree (pcm->dmodelsConv); + psFree (pcm->psf); + psFree (pcm->constraint); return; } -pmPCMData *pmPCMDataAlloc ( +pmPCMdata *pmPCMdataAlloc ( const psVector *params, const psVector *paramMask, pmSource *source) { - pmPCMData *pcm = (pmPCMData *) psAlloc(sizeof(pmPCMData)); - psMemSetDeallocator(pcm, (psFreeFunc) pmPCMDataFree); + pmPCMdata *pcm = (pmPCMdata *) psAlloc(sizeof(pmPCMdata)); + psMemSetDeallocator(pcm, (psFreeFunc) pmPCMdataFree); // Allocate storage images for raw model and derivative images - pcm->model = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - pcm->dmodels = psArrayAlloc (params->n); + pcm->modelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsFlux = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - pcm->dmodels->data[n] = NULL; + pcm->dmodelsFlux->data[n] = NULL; if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsFlux->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } // Allocate storage images for convolved model and derivative images - pcm->modelConv = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - pcm->dmodelsConv = psArrayAlloc (params->n); + pcm->modelConvFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsConvFlux = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - pcm->dmodelsConv->data[n] = NULL; + pcm->dmodelsConvFlux->data[n] = NULL; if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; } - pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - } + pcm->dmodelsConvFlux->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + } + + pcm->modelConv = NULL; + pcm->psf = NULL; + pcm->constraint = NULL; + pcm->nDOF = 0; return pcm; } -psImage *pmPCMDataSaveImage (pmPCMData *pcm) { - - psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32); +psImage *pmPCMdataSaveImage (pmPCMdata *pcm) { + + psImage *model = psImageCopy (NULL, pcm->modelConvFlux, PS_TYPE_F32); return model; } +psKernel *pmPCMkernelFromPSF (pmSource *source, int nPix) { + + assert (source); + assert (source->psfImage); // XXX build if needed? + + int x0 = source->peak->xf - source->psfImage->col0; + int y0 = source->peak->yf - source->psfImage->row0; + + // need to decide on the size: dynamically? statically? + psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); + + // XXX we should just re-construct a PSF at this location + // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER); + + // if the realized PSF for this object does not cover the full kernel, give up for now + if (x0 + psf->xMin < 0) goto escape; + if (x0 + psf->xMax >= source->psfImage->numCols) goto escape; + if (y0 + psf->yMin < 0) goto escape; + if (y0 + psf->yMax >= source->psfImage->numRows) goto escape; + + double sum = 0.0; + for (int j = psf->yMin; j <= psf->yMax; j++) { + for (int i = psf->xMin; i <= psf->xMax; i++) { + double value = source->psfImage->data.F32[y0 + j][x0 + i]; + psf->kernel[j][i] = value; + sum += value; + } + } + assert (sum > 0.0); + + // psf must be normalized (integral = 1.0) + for (int i = 0; i < psf->image->numRows; i++) { + for (int j = 0; j < psf->image->numCols; j++) { + psf->image->data.F32[i][j] /= sum; + } + } + + return psf; + +escape: + psFree (psf); + return NULL; +} + +pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize) { + + // make sure we savep a cached copy of the psf flux + pmSourceCachePSF (source, maskVal); + + // convert the cached cached psf model for this source to a psKernel + psKernel *psf = pmPCMkernelFromPSF (source, psfSize); + if (!psf) return NULL; + +# if (USE_DELTA_PSF) + psImageInit (psf->image, 0.0); + psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; +# endif + + // allocate the model + pmModel *modelConv = pmModelAlloc(modelType); + if (!modelConv) { + psFree (psf); + return NULL; + } + + // count the number of unmasked pixels: + int nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // XXX are we doing the right thing with the mask? + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[i][j])) { + continue; + } + nPix ++; + } + } + + psVector *params = modelConv->params; + + // create the minimization constraints + psMinConstraint *constraint = psMinConstraintAlloc(); + constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); + constraint->checkLimits = modelConv->modelLimits; + + // set parameter mask based on fitting mode + int nParams = 0; + switch (fitOptions->mode) { + case PM_SOURCE_FIT_NORM: + // NORM-only model fits only source normalization (Io) + nParams = 1; + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + break; + case PM_SOURCE_FIT_PSF: + // PSF model only fits x,y,Io + nParams = 3; + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0; + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + nParams = params->n - 1; + psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + break; + case PM_SOURCE_FIT_INDEX: + // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + if (params->n == 7) { + nParams = 1; + } else { + nParams = 2; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0; + } + break; + case PM_SOURCE_FIT_NO_INDEX: + // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params + psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + if (params->n == 7) { + nParams = params->n - 1; + } else { + nParams = params->n - 2; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1; + } + break; + default: + psAbort("invalid fitting mode"); + } + + // generate PCM data storage structure + pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source); + + pcm->modelConv = modelConv; + pcm->psf = psf; + pcm->constraint = constraint; + + if (nPix < nParams + 1) { + psTrace ("psModules.objects", 4, "insufficient valid pixels\n"); + pcm->modelConv->flags |= PM_MODEL_STATUS_BADARGS; + abort (); + // XXX This should not be an abort!! + } + pcm->nPix = nPix; + pcm->nDOF = nPix - nParams - 1; + + return pcm; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h (revision 28692) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h (revision 28692) @@ -0,0 +1,88 @@ +/* @file pmPCMdata.h + * structures and functions to support PSF-convolved model fitting + * + * @author EAM, IfA + * + * @version $Revision: 1.29 $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2010 Institute for Astronomy, University of Hawaii + */ + +# ifndef PM_PCM_DATA_H +# define PM_PCM_DATA_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +/** pmPCMdata : PSF Convolved Model data storage structure + * + * Structure to carry the data needed to generate a PSF-convolve model fit + * + */ +// +typedef struct { + psImage *modelFlux; + psArray *dmodelsFlux; + psImage *modelConvFlux; + psArray *dmodelsConvFlux; + + pmModel *modelConv; + psKernel *psf; + psMinConstraint *constraint; + int nPix; + int nDOF; +} pmPCMdata; + +// structures & functions to support psf-convolved model fitting + +// psf-convolved model fitting +bool psphotModelWithPSF_LMM ( + psMinimization *min, + psImage *covar, + psVector *params, + psMinConstraint *constraint, + pmSource *source, + const psKernel *psf, + psMinimizeLMChi2Func func); + +psF32 psphotModelWithPSF_SetABX( + psImage *alpha, + psVector *beta, + 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); + +pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize); + +psImage *pmPCMdataSaveImage (pmPCMdata *pcm); + +psF32 pmPCM_SetABX( + psImage *alpha, + psVector *beta, + const psVector *params, + const psVector *paramMask, + pmPCMdata *pcm, + const pmSource *source); + +bool pmPCM_MinimizeChisq ( + psMinimization *min, + psImage *covar, + psVector *params, + pmSource *source, + pmPCMdata *pcm); + +bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); + +bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); + + +/// @} +# endif /* PM_PCM_DATA_H */ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c (revision 28692) @@ -108,5 +108,5 @@ // This function calculates the psf and aperture magnitudes - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); // raw PSF mag, AP mag + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal); // raw PSF mag, AP mag if (!status || isnan(source->apMag)) { psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c (revision 28692) @@ -1,9 +1,45 @@ -# include "psphotInternal.h" -# define USE_DELTA_PSF 0 +/* @file pmSourceFitPCM.c + * structures and functions to support PSF-convolved model fitting + * + * @author EAM, IfA + * + * @version $Revision: 1.29 $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2010 Institute for Astronomy, University of Hawaii + */ -// 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_MIN_TOL = 0.1; -static psF32 PM_SOURCE_FIT_MODEL_MAX_TOL = 2.0; +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAMaskWeight.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPCMdata.h" + +# define FACILITY "psModules.objects" // input source has both modelPSF and modelEXT. on successful exit, we set the @@ -11,99 +47,26 @@ // convolved model image. -// XXX need to generalize this -- number of fitted parameters must be flexible based on the fitOptions - -pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { +bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // 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, psfSize); - if (!psf) return NULL; - -# if (USE_DELTA_PSF) - psImageInit (psf->image, 0.0); - psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; -# endif - - // generate copy of the model - // XXX we could modify the parameter values or even the model - // here based on the observed seeing (some lookup table...) - - // use the source moments, etc to guess basic model parameters - pmModel *modelConv = pmSourceModelGuess (source, modelType); - if (!modelConv) { - psFree (psf); - return NULL; - } - - // adjust the pixels based on the footprint - float radius = psphotSetRadiusEXT (readout, source, markVal); - if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false; - - // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: - psEllipseShape psfShape; - psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; - psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; - psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; - psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); - - // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: - psEllipseShape extShape; - extShape.sx = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2; - extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY]; - extShape.sy = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2; - psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0); - - // decrease the initial guess ellipse by psf_minor axis: - psEllipseAxes extAxesMod; - extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor))); - extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor))); - extAxesMod.theta = extAxes.theta; - - psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod); - modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2; - modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy; - modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2; - - // increase the initial guess central intensity by 2pi r^2: - modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor)); - - psVector *params = modelConv->params; - psVector *dparams = modelConv->dparams; - - // create the minimization constraints - psMinConstraint *constraint = psMinConstraintAlloc(); - constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); - constraint->checkLimits = modelConv->modelLimits; - - // set parameter mask based on fitting mode - // we fit a model without a floating sky term - int nParams = params->n - 1; - psVectorInit (constraint->paramMask, 0); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + psVector *params = pcm->modelConv->params; + psVector *dparams = pcm->modelConv->dparams; // force the floating parameters to fall within the contraint ranges for (int i = 0; i < params->n; i++) { - modelConv->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); - modelConv->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); + pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); + pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); } // set up the minimization process - psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_MIN_TOL, PM_SOURCE_FIT_MODEL_MAX_TOL); + psMinimization *myMin = psMinimizationAlloc (fitOptions->nIter, fitOptions->minTol, fitOptions->maxTol); psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - bool fitStatus = psphotModelWithPSF_LMM (myMin, covar, params, constraint, source, psf, modelConv->modelFunc); + bool fitStatus = pmPCM_MinimizeChisq (myMin, covar, params, source, pcm); for (int i = 0; i < dparams->n; i++) { - if (psTraceGetLevel("psphot") >= 4) { - fprintf (stderr, "%f ", params->data.F32[i]); - } - if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) + if ((pcm->constraint->paramMask != NULL) && pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); + psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]); } psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); @@ -118,12 +81,11 @@ // save the resulting chisq, nDOF, nIter - modelConv->chisq = myMin->value; - modelConv->nIter = myMin->iter; - - // XXX I actually need to count the number of unmasked pixels here - modelConv->nDOF = source->pixels->numCols*source->pixels->numRows - nParams; - - modelConv->flags |= PM_MODEL_STATUS_FITTED; - if (!fitStatus) modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; + pcm->modelConv->chisq = myMin->value; + pcm->modelConv->nIter = myMin->iter; + pcm->modelConv->nPix = pcm->nPix; + pcm->modelConv->nDOF = pcm->nDOF; + pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; + pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED; + if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; // models can go insane: reject these @@ -133,13 +95,50 @@ onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0); onPic &= (params->data.F32[PM_PAR_YPOS] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; + if (!onPic) pcm->modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed? - psFree(psf); psFree(myMin); psFree(covar); - psFree(constraint); - return modelConv; + return true; } + +bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { + + pcm->modelConv->modelGuess(pcm->modelConv, source); + + // generate copy of the model + // XXX we could modify the parameter values or even the model + // here based on the observed seeing (some lookup table...) + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape psfShape; + psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape extShape; + extShape.sx = pcm->modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2; + extShape.sxy = pcm->modelConv->params->data.F32[PM_PAR_SXY]; + extShape.sy = pcm->modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0); + + // decrease the initial guess ellipse by psf_minor axis: + psEllipseAxes extAxesMod; + extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor))); + extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor))); + extAxesMod.theta = extAxes.theta; + + psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod); + pcm->modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2; + pcm->modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy; + pcm->modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2; + + // increase the initial guess central intensity by 2pi r^2: + pcm->modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor)); + + return true; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28692) @@ -76,5 +76,5 @@ // XXX masked region should be (optionally) elliptical -bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal) +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -160,5 +160,5 @@ // measure the contribution of included pixels if (mode & PM_SOURCE_PHOT_DIFFSTATS) { - pmSourceMeasureDiffStats (source, maskVal); + pmSourceMeasureDiffStats (source, maskVal, markVal); } Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h (revision 28692) @@ -52,9 +52,9 @@ bool pmSourceMagnitudesInit (psMetadata *config); -bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal); +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal); bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal); bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor); -bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal); +bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal); Index: /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28691) +++ /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28692) @@ -149,4 +149,5 @@ #include #include +#include // The following headers are from random locations, here because they cross bounds Index: /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28692) @@ -168,5 +168,4 @@ psphotExtendedSourceFits.c \ psphotKernelFromPSF.c \ - psphotPSFConvModel.c \ psphotFitSet.c \ psphotSourceFreePixels.c \ Index: /branches/eam_branches/ipp-20100621/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28692) @@ -12,12 +12,4 @@ #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs - -// pmPCMData : PSF Convolved Model data storage structure -typedef struct { - psImage *model; - psArray *dmodels; - psImage *modelConv; - psArray *dmodelsConv; -} pmPCMData; // top-level psphot functions @@ -293,33 +285,4 @@ bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise); -// structures & functions to support psf-convolved model fitting - -// psf-convolved model fitting -bool psphotModelWithPSF_LMM ( - psMinimization *min, - psImage *covar, - psVector *params, - psMinConstraint *constraint, - pmSource *source, - const psKernel *psf, - psMinimizeLMChi2Func func); - -psF32 psphotModelWithPSF_SetABX( - psImage *alpha, - psVector *beta, - 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); - int psphotKapaOpen (void); bool psphotKapaClose (void); @@ -463,3 +426,6 @@ bool psphotStackObjectsUnifyPosition (psArray *objects); +bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal); +pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); + #endif Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotApResid.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotApResid.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotApResid.c (revision 28692) @@ -459,5 +459,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - bool status = pmSourceMagnitudes (source, psf, photMode, maskVal); + bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c (revision 28692) @@ -186,5 +186,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); if (status && isfinite(source->apMag)) Nap ++; Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28692) @@ -8,4 +8,5 @@ static int NfitDBL = 0; static int NfitEXT = 0; +static int NfitPCM = 0; bool psphotFitInit (int nThreads) { @@ -440,4 +441,8 @@ pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { + if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { + psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); + } + pmSourceFitOptions options = *fitOptions; @@ -456,8 +461,4 @@ } - if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { - psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); - } - // for sersic models, use a grid search to choose an index, then float the params there if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { @@ -476,5 +477,17 @@ } -pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { +pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { + psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); + } + + float radius = psphotSetRadiusEXT (readout, source, markVal); + + // XXX note that this changes the source moments that are published... + // recalculate the source moments using the larger extended-source moments radius + // at this stage, skip Gaussian windowing, and do not clip pixels by S/N + // this uses the footprint to judge both radius and aperture? + if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false; pmSourceFitOptions options = *fitOptions; @@ -487,18 +500,17 @@ // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); - // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelType); - if (!EXT) { + pmPCMdata *pcm = pmPCMinit (source, &options, modelType, maskVal, psfSize); + if (!pcm) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); return NULL; } - - if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { - psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); - } + // XXX check for nDOF too small + + // use the source moments, etc to guess basic model parameters + pmSourceModelGuessPCM (pcm, source, maskVal, markVal); // for sersic models, use a grid search to choose an index, then float the params there if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { - psphotFitSersicIndexPCM (source, EXT, fitOptions, maskVal); + psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal); } @@ -508,8 +520,8 @@ options.mode = PM_SOURCE_FIT_EXT; } - pmSourceFitPCM (source, EXT, &options, maskVal); + pmSourceFitPCM (source, PCM, &options, maskVal); // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); - return (EXT); + return (PCM); } @@ -520,5 +532,7 @@ // A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations -bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) { +bool psphotFitSersicIndex (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { + + pmModel *model = pcm->modelConv; assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); @@ -535,6 +549,7 @@ for (int i = 0; i < N_INDEX_GUESS; i++) { model->params->data.F32[PM_PAR_7] = indexGuess[i]; - model->modelGuess(model, source); - pmSourceFitModel (source, model, &options, maskVal); + pmSourceModelGuessPCM (pcm, source, maskVal, markVal); + + pmSourceFitPCM (pcm, source, &options, maskVal); chiSquare[i] = model->chisq; if (i == 0) { @@ -570,4 +585,6 @@ float xMin, chiSquare[N_INDEX_GUESS]; int iMin; + + // XXX we probably cannot be calling model->modelGuess() : this does not include the psf sigma for (int i = 0; i < N_INDEX_GUESS; i++) { Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28691) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28692) @@ -191,5 +191,5 @@ // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit @@ -342,5 +342,5 @@ // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit