Changeset 36085 for trunk/psModules/src/objects/pmPCMdata.c

Timestamp:

Aug 31, 2013, 5:55:16 AM (13 years ago)

Author:

eugene

Message:

merge changes from EAM branch (pmModel_CentralPixel functions for sersic-like model support; add interactive mode to PCM LMM fitting; harder PCM LMM damping (nu scaled by 3 not 2); add EXT_AND SKY and SHAPE fitting modes, include sky in INDEX-only fits; ifdefs for trace-only values; allow positions to go to 100,000; major re-work of central pixel & flux for sersic-like models; add pmPCMMakeModel function to generate the flux for an arbitrary model

Location:

trunk/psModules

Files:

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

trunk/psModules/src/objects/pmPCMdata.c

@@ -173 +173 @@
 }
 
+int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) {
+
+    // set parameter mask based on fitting mode
+    int nParams = 0;
+
+    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 = params->n - 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 = params->n;
+        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 (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:
+        // fits all but index (PAR7) including sky
+        psVectorInit (constraint->paramMask, 0);
+        if (params->n == 7) {
+            nParams = params->n;
+        } else {
+            nParams = params->n - 1;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+        }
+        break;
+      default:
+        psAbort("invalid fitting mode");
+    }
+    return nParams;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
 
@@ -219 +293 @@
     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");
-    }
+    int nParams = pmPCMsetParams (constraint, fitOptions->mode);
 
     if (nPix <  nParams + 1) {
@@ -341 +368 @@
     }
 
-    // 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) {
@@ -478 +457 @@
 }
 
+// construct a realization of the source model
+bool pmPCMMakeModel (pmSource *source, pmModel *model, 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);
+
+    // 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);
+    } 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;
+}
+