Changeset 37633

Timestamp:

Nov 19, 2014, 1:02:03 PM (12 years ago)

Author:

bills

Message:

add code to measure the Blakeslee bumpiness paramter

File:

• trunk/psphot/src/psphotGalaxyParams.c (modified) (8 diffs)

trunk/psphot/src/psphotGalaxyParams.c

@@ -10 +10 @@
 
     fprintf (stdout, "\n");
-    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Paramters ---");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Parameters ---");
 
     // select the appropriate recipe information
@@ -33 +33 @@
     for (int i = 0; i < num; i++) {
         if (!psphotGalaxyParamsReadout (config, view, filerule, i, recipe)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy paramter measurements for %s entry %d", filerule, i);
+            psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy parameter measurements for %s entry %d", filerule, i);
             return false;
         }
@@ -256 +256 @@
     // don't try and use bad models
     pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_NONCONVERGE;
     badModel |= PM_MODEL_STATUS_BADARGS;
     badModel |= PM_MODEL_STATUS_OFFIMAGE;
@@ -337 +338 @@
             // do we need this?
             source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-            if (source->modelEXT->isPCM) {
-                pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
-            } else {
-                pmSourceCacheModel (source, maskVal);
-            }
-        }
+#ifdef notdef
+#endif
+        }
+
 
         psF32 *PAR = sersicModel->params->data.F32;
         psF32 Xo = PAR[PM_PAR_XPOS];
         psF32 Yo = PAR[PM_PAR_YPOS];
-        
+        psF32 Io = PAR[PM_PAR_I0];
         // For now: set the half light radius to the major axis of the model
         psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff);
@@ -363 +362 @@
         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius);
 
+        // RedefinePixels frees modelFlux
+        if (source->modelEXT->isPCM) {
+            pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
+        } else {
+            pmSourceCacheModel (source, maskVal);
+        }
 
 #ifdef DUMP_IMAGES
@@ -485 +490 @@
         vPix = source->pixels->data.F32;
         vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+        psF32 **vVar = source->variance->data.F32;
+        psF32 **vModel = source->modelFlux->data.F32;
         // Accumulate Rt and Ra
         psF32 sumRt = 0;
         psF32 sumRa = 0;
         int numPixels = 0;
+        psF32 sumB = 0;
+        psF32 sumModel = 0;
+        int numB = 0;
         for (psS32 row = 0; row < source->pixels->numRows ; row++) {
             psF32 yDiff = row - yCM;
@@ -506 +516 @@
                 if (r2 > R2) continue;
 
+                // Measurement of bumpiness parameter excludes the pixels
+                // within 2 pixels of the center
+                if (r2 > 4) {
+                    psF32 variance = vVar[row][col];
+                    if (isfinite(variance)) {
+                        psF32 modelPixel = vModel[row][col];
+                        if (!isfinite(modelPixel)) {
+                            // XXX: Can this happen in a good model?
+                            fprintf(stderr, "non-finite modelPixel %4d (%4d, %4d)\n", source->id, col, row);
+                        }
+                        numB++;
+                        sumB += PS_SQR(pixel) - variance;
+                        sumModel += Io * modelPixel;
+                    }
+                }
+
                 // x coordinate of mirror pixel
                 int xFlip = xCM - xDiff;
@@ -519 +545 @@
                 ++numPixels;
             }
+        }
+
+        if (numB > 0) {
+            source->extpars->gbumpy = sqrt(sumB/numB)  / (sumModel / numB);
         }