Changeset 28692

Timestamp:

Jul 20, 2010, 7:21:53 PM (16 years ago)

Author:

eugene

Message:

update PCM fitting (not quite ready)

Location:

branches/eam_branches/ipp-20100621

Files:

• psModules/src/objects/Makefile.am (modified) (2 diffs)
• psModules/src/objects/pmPCM_MinimizeChisq.c (modified) (13 diffs)
• psModules/src/objects/pmPCMdata.c (modified) (1 diff)
• psModules/src/objects/pmPCMdata.h (added)
• psModules/src/objects/pmPSFtryFitPSF.c (modified) (1 diff)
• psModules/src/objects/pmSourceFitPCM.c (modified) (4 diffs)
• psModules/src/objects/pmSourcePhotometry.c (modified) (2 diffs)
• psModules/src/objects/pmSourcePhotometry.h (modified) (1 diff)
• psModules/src/psmodules.h (modified) (1 diff)
• psphot/src/Makefile.am (modified) (1 diff)
• psphot/src/psphot.h (modified) (3 diffs)
• psphot/src/psphotApResid.c (modified) (1 diff)
• psphot/src/psphotMagnitudes.c (modified) (1 diff)
• psphot/src/psphotSourceFits.c (modified) (9 diffs)
• psphot/src/psphotSourceSize.c (modified) (2 diffs)

branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am

@@ -61 +61 @@
         pmPSFtryFitPSF.c \
         pmPSFtryMetric.c \
+        pmPCMdata.c \
+        pmPCM_MinimizeChisq.c \
+        pmSourceFitPCM.c \
         pmTrend2D.c \
         pmGrowthCurveGenerate.c \
@@ -105 +108 @@
         pmPSF_IO.h \
         pmPSFtry.h \
+        pmPCMdata.h \
         pmTrend2D.h \
         pmGrowthCurve.h \

branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -1 @@
-# 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 <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <strings.h>
+#include <pslib.h>
+#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 +86 @@
     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 +93 @@
     }
     // 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 +116 @@
 
         // 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 +140 @@
         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 +149 @@
 
         /* 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 +159 @@
             // 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 +202 @@
     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 +243 @@
     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 +278 @@
             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 +291 @@
 
     // 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 +332 @@
             }
 
-            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 +344 @@
             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++;

branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c

@@ -1 @@
-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 <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <strings.h>
+#include <pslib.h>
+#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;
+}

branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c

@@ -108 +108 @@
 
         // 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

branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c

@@ -1 @@
-# 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 <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <strings.h>
+#include <pslib.h>
+#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 +47 @@
 // 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 +81 @@
 
     // 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 +95 @@
     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;
+}

branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c

@@ -76 +76 @@
 
 // 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 +160 @@
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
-        pmSourceMeasureDiffStats (source, maskVal);
+        pmSourceMeasureDiffStats (source, maskVal, markVal);
     }
 

branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.h

@@ -52 +52 @@
 
 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);

branches/eam_branches/ipp-20100621/psModules/src/psmodules.h

@@ -149 +149 @@
 #include <pmSourceMatch.h>
 #include <pmDetEff.h>
+#include <pmPCMdata.h>
 
 // The following headers are from random locations, here because they cross bounds

branches/eam_branches/ipp-20100621/psphot/src/Makefile.am

@@ -168 +168 @@
         psphotExtendedSourceFits.c     \
         psphotKernelFromPSF.c          \
-        psphotPSFConvModel.c           \
         psphotFitSet.c                 \
         psphotSourceFreePixels.c       \

branches/eam_branches/ipp-20100621/psphot/src/psphot.h

@@ -12 +12 @@
 
 #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 +285 @@
 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 +426 @@
 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

branches/eam_branches/ipp-20100621/psphot/src/psphotApResid.c

@@ -459 +459 @@
         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

branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c

@@ -186 +186 @@
         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 ++;
 

branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c

@@ -8 +8 @@
 static int NfitDBL = 0;
 static int NfitEXT = 0;
+static int NfitPCM = 0;
 
 bool psphotFitInit (int nThreads) {
@@ -440 +441 @@
 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 +461 @@
     }
 
-    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 +477 @@
 }
 
-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 +500 @@
     // 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 +520 @@
         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 +532 @@
 // 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 +549 @@
     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 +585 @@
     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++) {

branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c

@@ -191 +191 @@
 
         // 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 +342 @@
 
         // 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