Changeset 33415 for branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

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

branches/meh_branches/ppstack_test/psModules/src/objects

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c
+pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V3.v1.c
+pmSourceIO_CMF_PS1_V1.v1.c
+pmSourceIO_CMF_PS1_V2.v1.c
+

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c
+pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V3.v1.c
+pmSourceIO_CMF_PS1_V1.v1.c
+pmSourceIO_CMF_PS1_V2.v1.c
+

branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c

@@ -41 +41 @@
 #include "pmPCMdata.h"
 
+# define USE_DELTA_PSF 0
+# define USE_1D_GAUSS 1
+
 static void pmPCMdataFree (pmPCMdata *pcm) {
 
@@ -54 +57 @@
     psFree (pcm->psfFFT);
     psFree (pcm->constraint);
+    psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -88 +92 @@
     pcm->constraint = NULL;
     pcm->nDOF = 0;
+
+    // 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->sigma = 1.0; // this should be set to something sensible when the psf is known
 
     return pcm;
@@ -257 +267 @@
     pcm->nDOF = nPix - nParams;
 
+# if (USE_1D_GAUSS)
+    pmModel *modelPSF = source->modelPSF;
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+    axes = psEllipseShapeToAxes (shape, 20.0);
+    
+    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[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);
+# else
+    pcm->smdata = NULL;
+# endif
+
     return pcm;
 }
@@ -364 +398 @@
             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);
     }
 
     return true;
 }
+
+// construct a realization of the source model
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    // 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");
+    
+        psEllipseShape shape;
+        psEllipseAxes axes;
+
+        shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+        shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+        shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+        axes = psEllipseShapeToAxes (shape, 20.0);
+    
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[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;
+}
+

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c
+pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V3.v1.c
+pmSourceIO_CMF_PS1_V1.v1.c
+pmSourceIO_CMF_PS1_V2.v1.c
+