Changeset 37066 for branches/eam_branches/ipp-ops-20130712/psModules/src/objects/pmPCMdata.c

Timestamp:

Jul 17, 2014, 12:30:45 PM (12 years ago)

Author:

eugene

Message:

merge changes (from past YEAR) into trunk

Location:

branches/eam_branches/ipp-ops-20130712/psModules

Files:

• . (modified) (1 prop)
• src/objects (modified) (1 prop)
• src/objects/pmPCMdata.c (modified) (14 diffs)

branches/eam_branches/ipp-ops-20130712/psModules/src/objects

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c
+pmSourceIO_CMF_PS1_SV3.c
 pmSourceIO_CMF_PS1_DV1.c
 pmSourceIO_CMF_PS1_DV2.c
 pmSourceIO_CMF_PS1_DV3.c
-
+pmSourceIO_CMF_PS1_DV4.c

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c
+pmSourceIO_CMF_PS1_SV3.c
 pmSourceIO_CMF_PS1_DV1.c
 pmSourceIO_CMF_PS1_DV2.c
 pmSourceIO_CMF_PS1_DV3.c
-
+pmSourceIO_CMF_PS1_DV4.c

branches/eam_branches/ipp-ops-20130712/psModules/src/objects/pmPCMdata.c

@@ -31 +31 @@
 #include "pmMoments.h"
 #include "pmModelFuncs.h"
+#include "pmModelClass.h"
 #include "pmModel.h"
 #include "pmModelUtils.h"
-#include "pmModelClass.h"
 #include "pmSourceMasks.h"
 #include "pmSourceExtendedPars.h"
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -43 +44 @@
 
 # define USE_DELTA_PSF 0
-# define USE_1D_GAUSS 1
 
 static void pmPCMdataFree (pmPCMdata *pcm) {
@@ -58 +58 @@
     psFree (pcm->psfFFT);
     psFree (pcm->constraint);
+
     psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -88 +90 @@
     }
 
+    pcm->smdata = NULL;
+    pcm->smdata2d = NULL;
+
     pcm->modelConv = NULL;
     pcm->psf = NULL;
@@ -94 +99 @@
     pcm->nDOF = 0;
 
+    pcm->poissonErrors = true;
+
     // full convolution with the PSF is expensive.  if we have to save time, we can do a 1D
     // convolution with a Gaussian approximation to the kernel
     pcm->use1Dgauss = false;
-    pcm->nsigma = 3.0; 
+    pcm->nsigma = NAN; // this is set to something defined by the user
     pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
 
@@ -173 +180 @@
 }
 
+int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) {
+
+    // set parameter mask based on fitting mode
+    int nParams = 0;
+    int nParAll = constraint->paramMask->n;
+
+    switch (mode) {
+      case PM_SOURCE_FIT_NORM:
+        // 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:
+        // fits only 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:
+        // fits all params except sky
+        nParams = nParAll - 1;
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+        break;
+
+      case PM_SOURCE_FIT_EXT_AND_SKY:
+        // fits all params including sky
+        nParams = nParAll;
+        psVectorInit (constraint->paramMask, 0);
+        break;
+
+      case PM_SOURCE_FIT_SHAPE:
+        // fits shape (Sxx, Sxy, Syy) and Io
+        nParams = 5;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXX] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXY] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SYY] = 0;
+        break;
+
+      case PM_SOURCE_FIT_INDEX:
+        // fits only 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 (nParAll == 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:
+        // fits all but index (PAR7) including sky
+        psVectorInit (constraint->paramMask, 0);
+        if (nParAll == 7) {
+            nParams = nParAll;
+        } else {
+            nParams = nParAll - 1;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+        }
+        break;
+      default:
+        psAbort("invalid fitting mode");
+    }
+    return nParams;
+}
+
+static int modelType_GAUSS = -1;
+static int modelType_PS1_V1 = -1;
+
+// generate a Gaussian smoothing kernel for supplied sigma.  sigma here does not need to match
+// that used to allocate the structure, but it is recommended
+bool psImageSmoothCacheKernel_PS1_V1 (psImageSmoothCacheData *smdata, float sigma, float kappa) {
+    // check for NULL structure elements?
+
+    int size = smdata->Nrange;
+
+    psFree (smdata->kernel);
+    smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32);
+
+    double sum = 0.0;                   // Sum of Gaussian, for normalization
+    double factor = 1.0 / (sigma * M_SQRT2);    // Multiplier for i -> z
+
+    // PS1_V1 is a power-law with fitted linear term:
+    // 1 / (1 + kappa z + z^1.666)  where z = (r/sigma)^2
+
+    // generate the kernel (not normalized)
+    for (int i = -size, j = 0; i <= size; i++, j++) {
+        float z = PS_SQR(i * factor);
+        sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666));
+    }
+
+    // renormalize kernel to integral of 1.0
+    for (int i = 0; i < 2 * size + 1; i++) {
+        smdata->kernel->data.F32[i] /= sum;
+    }
+
+    return true;
+}
+
+psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
+
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseAxes axes;
+    bool useReff = modelPSF->class->useReff;
+    psF32 *PAR = modelPSF->params->data.F32;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+        float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        *sigma = 0.5 * (axes.major + axes.minor);
+        *kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (flux, *sigma, nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+        psImageSmoothCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
+psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
+
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseAxes axes;
+    bool useReff = modelPSF->class->useReff;
+    psF32 *PAR = modelPSF->params->data.F32;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+        float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        *sigma = 0.5 * (axes.major + axes.minor);
+        *kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+        psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
 
-    // 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 = pmPCMkernelFromPSF (source, psfSize);
-    if (!psf) {
-        // NOTE: this only happens if the source is too close to an edge
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        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
+    modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS");
+    modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1");
 
     // count the number of unmasked pixels:
@@ -217 +393 @@
     psMinConstraint *constraint = psMinConstraintAlloc();
     constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
-    constraint->checkLimits = model->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");
-    }
+    constraint->checkLimits = model->class->modelLimits;
+
+    int nParams = pmPCMsetParams (constraint, fitOptions->mode);
 
     if (nPix <  nParams + 1) {
         psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
-        psFree (psf);
         psFree (constraint);
         model->flags |= PM_MODEL_STATUS_BADARGS;
@@ -278 +406 @@
     // generate PCM data storage structure
     pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
-
-    pcm->psf = psf;
     pcm->modelConv = psMemIncrRefCounter(model);
     pcm->constraint = constraint;
+
+    pcm->poissonErrors = fitOptions->poissonErrors;
+    pcm->nsigma = fitOptions->nsigma;
 
     pcm->nPix = nPix;
@@ -288 +417 @@
 
 # if (USE_1D_GAUSS)
-    pmModel *modelPSF = source->modelPSF;
-    psAssert (modelPSF, "psf model must be defined");
-    
-    psEllipseAxes axes;
-    bool useReff = pmModelUseReff (modelPSF->type);
-    psF32 *PAR = modelPSF->params->data.F32;
-    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
-    
-    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
-    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
 
     pcm->use1Dgauss = true;
-    pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-    pcm->nsigma = 2.0;
-
-    pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+    if (USE_1D_CACHE) {
+        pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    } else {
+        pcm->smdata2d = psImageSmooth2dCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    }
+
 # else
+    // 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 = pmPCMkernelFromPSF (source, psfSize);
+    if (!psf) {
+        // NOTE: this only happens if the source is too close to an edge
+        model->flags |= PM_MODEL_STATUS_BADARGS;
+        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
+    pcm->psf = psf;
     pcm->smdata = NULL;
 # endif
@@ -341 +478 @@
     }
 
-    // if we changed the fit mode, we need to update nDOF
-    int nParams = 0;
-    // set parameter mask based on fitting mode
-    switch (fitOptions->mode) {
-      case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = 1;
-        psVectorInit (pcm->constraint->paramMask, 1);
-        pcm->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 (pcm->constraint->paramMask, 1);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
-        pcm->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 = model->params->n - 1;
-        psVectorInit (pcm->constraint->paramMask, 0);
-        pcm->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 (pcm->constraint->paramMask, 1);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        if (model->params->n == 7) {
-            nParams = 1;
-        } else {
-            nParams = 2;
-            pcm->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 (pcm->constraint->paramMask, 0);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        if (model->params->n == 7) {
-            nParams = model->params->n - 1;
-        } else {
-            nParams = model->params->n - 2;
-            pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
-        }
-        break;
-      default:
-        psAbort("invalid fitting mode");
-    }
+    int nParams = pmPCMsetParams (pcm->constraint, fitOptions->mode);
 
     if (pcm->nPix <  nParams + 1) {
@@ -415 +504 @@
             pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
         }
-        psFree(pcm->smdata);
-        pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+        // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
+        if (USE_1D_CACHE) {
+            psFree(pcm->smdata);
+            pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+            pmModel *modelPSF = source->modelPSF;
+            if (modelPSF->type == modelType_GAUSS) {
+                psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
+            }
+            if (modelPSF->type == modelType_PS1_V1) {
+                psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
+            }
+        } else {
+            psFree(pcm->smdata2d);
+            pcm->smdata2d = psImageSmooth2dCacheAlloc (pcm->nsigma);
+
+            pmModel *modelPSF = source->modelPSF;
+            if (modelPSF->type == modelType_GAUSS) {
+                // psImageSmooth2dCacheKernel_Gauss (pcm->smdata2d, pcm->sigma);
+            }
+            if (modelPSF->type == modelType_PS1_V1) {
+                psImageSmooth2dCacheKernel_PS1_V1 (pcm->smdata2d, pcm->sigma, pcm->kappa);
+            }
+        }
     }
 
@@ -423 +535 @@
 
 // construct a realization of the source model
-bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -440 +552 @@
     // convolve the model image with the PSF
     if (USE_1D_GAUSS) {
-        // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
-        // * the model flux is not masked
-        // * threading takes place above this level
         
-        // define the Gauss parameters from the psf
-        pmModel *modelPSF = source->modelPSF;
-        psAssert (modelPSF, "psf model must be defined");
-    
-        psEllipseAxes axes;
-        bool useReff = pmModelUseReff (modelPSF->type);
-        psF32 *PAR = modelPSF->params->data.F32;
-        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
-    
-        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
-        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-
-        float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-        float nsigma = 2.0;
-
-        psImageSmooth (source->modelFlux, sigma, nsigma);
+        float sigma = NAN;
+        float kappa = NAN;
+
+        if (USE_1D_CACHE) {
+            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+            psImageSmoothCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        } else {
+            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        }
+        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
     } else {
         // make sure we save a cached copy of the psf flux
@@ -478 +584 @@
 }
 
+// construct a realization of the source model
+bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // modelFlux always has unity normalization (I0 = 1.0)
+    // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | PM_MODEL_OP_NORM, maskVal);
+
+    // convolve the model image with the PSF
+    if (USE_1D_GAUSS) {
+
+        float sigma = NAN;
+        float kappa = NAN;
+
+        if (USE_1D_CACHE) {
+            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+            psImageSmoothCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        } else {
+            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        }
+        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
+    } else {
+        // 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 = pmPCMkernelFromPSF (source, psfSize);
+        if (!psf) {
+            // NOTE: this only happens if the source is too close to an edge
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return NULL;
+        }
+
+        // XXX not sure if I can place the output on top of the input
+        psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf);
+    }
+    return true;
+}

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c
+pmSourceIO_CMF_PS1_SV3.c
 pmSourceIO_CMF_PS1_DV1.c
 pmSourceIO_CMF_PS1_DV2.c
 pmSourceIO_CMF_PS1_DV3.c
-
+pmSourceIO_CMF_PS1_DV4.c