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

Timestamp:

Nov 26, 2013, 5:51:56 AM (13 years ago)

Author:

eugene

Message:

fixing the non-poisson chisq for pcm fits; normalize sky if we as for a normalized model

Location:

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

Files:

• pmModel.c (modified) (1 diff)
• pmPCM_MinimizeChisq.c (modified) (1 diff)
• pmPCMdata.c (modified) (1 diff)
• pmSourceFitPCM.c (modified) (2 diffs)
• pmSourcePhotometry.c (modified) (1 diff)

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

@@ -217 +217 @@
     // the options allow us to modify various aspects of the model
     if (mode & PM_MODEL_OP_NORM) {
+        // if we are including the sky, renormalizing should force use to normalized down the sky flux
+        params->data.F32[PM_PAR_SKY] /= params->data.F32[PM_PAR_I0];
         params->data.F32[PM_PAR_I0] = 1.0;
     }
     if (!(mode & PM_MODEL_OP_SKY)) {
         params->data.F32[PM_PAR_SKY] = 0.0;
-    }
+    } 
     if (mode & PM_MODEL_OP_CENTER) {
         params->data.F32[PM_PAR_XPOS] = image->col0 + 0.5*image->numCols;

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

@@ -135 +135 @@
         }
 
-        if (1) {
+        if (0) {
             fprintf (stderr, "%d : ", min->iter);
             for (int ti = 0; ti < params->n; ti++) {

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

@@ -597 +597 @@
 
     // modelFlux always has unity normalization (I0 = 1.0)
-    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | PM_MODEL_OP_NORM, maskVal);
 
     // convolve the model image with the PSF

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

@@ -51 +51 @@
 # define TIMING 0
 
+bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal);
+
 bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     
@@ -113 +115 @@
     } else {
         // xxx this is wrong because it does not convolve with the psf
-        pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
+        pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+        pmSourceChisqModelFlux (source, pcm->modelConv, maskVal);
     }
     if (TIMING) { t4 = psTimerMark ("pmSourceFitPCM"); }

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

@@ -902 +902 @@
 }
 
+bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+
+    float dC = 0.0;
+    int Npix = 0;
+
+    psVector *params = model->params;
+    psImage  *image = source->pixels;
+    psImage  *modelFlux = source->modelFlux;
+    psImage  *mask = source->maskObj;
+    psImage  *variance = source->variance;
+
+    float Io = params->data.F32[PM_PAR_I0];
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+
+            // skip pixels which are masked
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+
+            if (variance->data.F32[iy][ix] <= 0) continue;
+
+            dC += PS_SQR (image->data.F32[iy][ix] - Io*modelFlux->data.F32[iy][ix]) / variance->data.F32[iy][ix];
+            Npix ++;
+        }
+    }
+    model->nPix = Npix;
+    model->nDOF = Npix - model->nPar;
+    model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
+
+    return (true);
+}
+
 double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
 {