Changeset 36021 for branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c

Timestamp:

Aug 25, 2013, 12:08:56 PM (13 years ago)

Author:

eugene

Message:

convert SERSIC to use pmModelCP_SersicSubpix; clean up DEV and EXP versions; add pmPCMMakeModel, pmSersicKappa, pmSersicNorm functions

File:

• branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c (modified) (2 diffs)

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

@@ -261 +261 @@
         // 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;
+        // constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         if (params->n == 7) {
+            nParams = params->n;
+        } else {
             nParams = params->n - 1;
-        } else {
-            nParams = params->n - 2;
             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
         }
@@ -488 +488 @@
 }
 
+// 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;
+}
+