Changeset 36229 for branches/eam_branches/ipp-20130904/psModules

Timestamp:

Oct 22, 2013, 10:30:55 AM (13 years ago)

Author:

eugene

Message:

swap to PS1_V1 smoothing for the pcm analysis

Location:

branches/eam_branches/ipp-20130904/psModules/src/objects

Files:

• pmPCM_MinimizeChisq.c (modified) (3 diffs)
• pmPCMdata.c (modified) (8 diffs)
• pmPCMdata.h (modified) (1 diff)

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

@@ -477 +477 @@
         // * threading takes place above this level
         pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
-        psImageSmooth_PreAlloc_F32 (pcm->modelConvFlux, pcm->smdata);
-        // psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma);
+        psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata);
     } else {
         psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
@@ -497 +496 @@
             // * threading takes place above this level
             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
-            psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
-            // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
+            psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
         } else {
             psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
@@ -518 +516 @@
             // * threading takes place above this level
             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
-            psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
-            // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
+            psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
         } else {
             psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);

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

@@ -88 +88 @@
     }
 
+    pcm->smdata = NULL;
     pcm->modelConv = NULL;
     pcm->psf = NULL;
@@ -253 +254 @@
     // check for NULL structure elements?
 
-    int size = smdata->Nrange + 1;
-
+    int size = smdata->Nrange;
+
+    psFree (smdata->kernel);
     smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32);
 
@@ -265 +267 @@
     // generate the kernel (not normalized)
     for (int i = -size, j = 0; i <= size; i++, j++) {
-        float z = i / sigma / M_SQRT2;
+        float z = PS_SQR(i * factor);
         sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666));
     }
@@ -334 +336 @@
     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;
+    
+    // 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);
 
-    float v1sigma = 0.5 * (axes.major + axes.minor);
-    psImageSmoothCacheKernel_Gauss (pcm->smdata, v1sigma, PAR[PM_PAR_7]);
-
+    // 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);
 # else
     // make sure we save a cached copy of the psf flux
@@ -363 +372 @@
     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
 
@@ -428 +435 @@
             pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
         }
+
+        // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
         psFree(pcm->smdata);
-        pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
-
-    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
-    pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
-
-    float v1sigma = 0.5 * (axes.major + axes.minor);
-    psImageSmoothCacheKernel_Gauss (pcm->smdata, v1sigma, PAR[PM_PAR_7]);
+        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);
     }
 
@@ -472 +481 @@
         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);
+        // 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
-        pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
-        
-        float v1sigma = 0.5 * (axes.major + axes.minor);
-        psImageSmoothCacheKernel_Gauss (pcm->smdata, v1sigma, PAR[PM_PAR_7]);
-
-        psImageSmoothCache_F32 (source->modelFlux, pcm->smdata);
+        psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (source->modelFlux, sigma, nsigma);
+
+        psImageSmoothCacheKernel_PS1_V1 (smdata, sigma, kappa);
+        // XXX psImageSmoothCacheKernel_GAUSS (smdata, sigma);
+
+        psImageSmoothCache_F32 (source->modelFlux, smdata);
+
+        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
     } else {
         // make sure we save a cached copy of the psf flux
@@ -532 +545 @@
         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 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);
+
+        psImageSmoothCache_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

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

@@ -38 +38 @@
     bool use1Dgauss;
     float sigma;
+    float kappa;
     float nsigma;
 
-    psImageSmooth_PreAlloc_Data *smdata;
+    psImageSmoothCacheData *smdata;
 } pmPCMdata;