Changeset 36235 for branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c

Timestamp:

Oct 23, 2013, 10:31:13 AM (13 years ago)

Author:

eugene

Message:

allow PS1_V1 or GAUSS smoothing depending on the psf model

File:

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

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

@@ -249 +249 @@
 }
 
+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
@@ -279 +282 @@
 }
 
+psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, psImage *flux, pmModel *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 nsigma = 2.0;
+    *sigma = NAN;
+    *kappa = NAN;
+
+    if (modelPSF->type == modelType_GAUSS) {
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        *sigma = 0.5 * (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;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
+
+    modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS");
+    modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1");
 
     // count the number of unmasked pixels:
@@ -328 +371 @@
 
 # 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);
-    
+
     pcm->use1Dgauss = true;
-    pcm->nsigma = 2.0;
-    
-    // XXX use these lines for a GAUSS kernel:
-    // float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
-    // float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-    // pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-
-    // XXX use these lines for a PS1_V1 kernel
-    pcm->sigma = 0.5 * (axes.major + axes.minor);
-    pcm->kappa = PAR[PM_PAR_7];
-
-    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
-    pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
-
-    // XXX use this for GAUSS kernel:
-    // psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
-
-    // XXX use these lines for a PS1_V1 kernel
-    psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
+    pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, source->pixels, source->modelPSF);
+
 # else
     // make sure we save a cached copy of the psf flux
@@ -440 +459 @@
         pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
 
-        // XXX use this for GAUSS kernel:
-        // psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
-
-        // XXX use these lines for a PS1_V1 kernel
-        psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
+        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);
+        }
     }
 
@@ -468 +489 @@
     // 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 nsigma = 2.0;
-
-        // XXX use these for a GAUSS kernel
-        // 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 sigma = 0.5 * (axes.major + axes.minor);
-        float kappa = PAR[PM_PAR_7];
-
-        // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
-        psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (source->modelFlux, sigma, nsigma);
-
-        psImageSmoothCacheKernel_PS1_V1 (smdata, sigma, kappa);
-        // XXX psImageSmoothCacheKernel_GAUSS (smdata, sigma);
+        float sigma = NAN;
+        float kappa = NAN;
+        psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, source->modelFlux, source->modelPSF);
 
         psImageSmoothCache_F32 (source->modelFlux, smdata);
-        psFree(smdata);
+        psFree (smdata);
 
         // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
@@ -533 +530 @@
     // 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 nsigma = 2.0;
-
-        // 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 sigma = 0.5 * (axes.major + axes.minor);
-        float kappa = PAR[PM_PAR_7];
-
-        // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
-        psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (source->modelFlux, sigma, nsigma);
-
-        psImageSmoothCacheKernel_PS1_V1 (smdata, sigma, kappa);
-        // XXX psImageSmoothCacheKernel_GAUSS (smdata, sigma);
+
+        float sigma = NAN;
+        float kappa = NAN;
+        psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, source->modelFlux, source->modelPSF);
 
         psImageSmoothCache_F32 (source->modelFlux, smdata);
@@ -582 +556 @@
     return true;
 }
-