Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/imcombine/pmPSFEnvelope.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmPSFEnvelope.c (modified) (10 diffs)

branches/simtest_nebulous_branches/psModules/src/imcombine/pmPSFEnvelope.c

@@ -34 +34 @@
 
 // #define TESTING                         // Enable test output
+// #define PEAK_NORM                       // Normalise peaks?
 #define PEAK_FLUX 1.0e4                 // Peak flux for each source
 #define SKY_VALUE 0.0e0                 // Sky value for fake image
@@ -64 +65 @@
                      int radius,        // Radius of each PSF
                      const char *modelName,// Name of PSF model to use
-                     int xOrder, int yOrder // Order for PSF variation fit
+                     int xOrder, int yOrder, // Order for PSF variation fit
+                     psImageMaskType maskVal
                      )
 {
@@ -122 +124 @@
             continue;
         }
+
+        if (psTraceGetLevel("psModules.imcombine") >= 1) {
+            psString string = NULL;     // String with values
+            psStringAppend(&string, "PSF %d: ", i);
+            float x = numCols / 2.0, y = numRows / 2.0; // Coordinates of interest
+            for (int j = 4; j < psf->params->n; j++) {
+                pmTrend2D *trend = psf->params->data[j]; // Trend of interest
+                double val = pmTrend2DEval(trend, x, y);
+                double err;
+                switch (trend->mode) {
+                  case PM_TREND_POLY_ORD:
+                  case PM_TREND_POLY_CHEB:
+                    err = NAN;
+                    break;
+                  case PM_TREND_MAP:
+                    err = psImageUnbinPixel(x, y, trend->map->error, trend->map->binning);
+                    break;
+                  default:
+                    psAbort("Unsupported mode: %x", trend->mode);
+                }
+                psStringAppend(&string, "%lf %lf   ", val, err);
+            }
+            psTrace("psModules.imcombine", 1, "%s\n", string);
+            psFree(string);
+        }
+
+        // Test PSF
+        {
+            bool goodPSF = false;       // Is there a PSF that we can use?
+            int xNum = PS_MAX(psf->trendNx, 1), yNum = PS_MAX(psf->trendNy, 1); // Number of positions to check
+            for (int j = 0; j < yNum && !goodPSF; j++) {
+                float y = ((float)j + 0.5) / (float)yNum * numRows; // Position on image
+                for (int i = 0; i < xNum && !goodPSF; i++) {
+                    float x = ((float)i + 0.5) / (float)xNum * numCols; // Position on image
+                    pmModelClassSetLimits(PM_MODEL_LIMITS_IGNORE);
+                    pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Test model
+                    if (!model) {
+                        continue;
+                    }
+                    model->modelSetLimits(PM_MODEL_LIMITS_MODERATE);
+                    bool limits = true; // Model within limits?
+                    for (int j = 0; j < model->params->n && limits; j++) {
+                        if (!model->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) ||
+                            !model->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) {
+                            limits = false;
+                        }
+                    }
+                    psFree(model);
+                    if (limits) {
+                        goodPSF = true;
+                    }
+                }
+            }
+            if (!goodPSF) {
+                psWarning("PSF %d is completely bad --- not including in envelope calculation.", i);
+                continue;
+            }
+        }
+
         pmResiduals *resid = psf->residuals;// PSF residuals
         psf->residuals = NULL;
-        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf,
-                                      NAN, radius, true, true)) {
+        pmModelClassSetLimits(PM_MODEL_LIMITS_MODERATE);
+        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
+                                      NAN, radius, true, false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
             psFree(envelope);
@@ -144 +206 @@
             float y = source->peak->yf + yOffset->data.S32[j]; // y coordinate of source
 
-            double flux = fakeRO->image->data.F32[(int)y][(int)x];
+#ifdef PEAK_NORM
+            // Perhaps I'm being paranoid, but specify a range to check
+            int uMax = PS_MIN(x + radius, numCols - 1), uMin = PS_MAX(x - radius, 0);
+            int vMax = PS_MIN(y + radius, numRows - 1), vMin = PS_MAX(y - radius, 0);
+
+            double flux = -INFINITY;    // Peak flux
+            for (int v = vMin; v <= vMax; v++) {
+                for (int u = uMin; u <= uMax; u++) {
+                    if (fakeRO->image->data.F32[v][u] > flux) {
+                        flux = fakeRO->image->data.F32[v][u];
+                    }
+                }
+            }
             if (!isfinite(flux) || flux < 0) {
                 continue;
             }
             float norm = PEAK_FLUX / flux; // Normalisation for source
+#endif
             psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region
             psImage *subImage = psImageSubset(fakeRO->image, region); // Subimage of fake PSF
             psImage *subEnv = psImageSubset(envelope, region); // Subimage of envelope
+#ifdef PEAK_NORM
             psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32));
+#endif
             psBinaryOp(subEnv, subEnv, "MAX", subImage);
             psFree(subImage);
@@ -163 +240 @@
             }
             float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source
+            psFree(model);
             if (srcRadius == 0) {
                 continue;
@@ -298 +376 @@
         }
 
-        // measure the source moments: tophat windowing, no pixel S/N cutoff
-        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
+        // measure the source moments: tophat windowing, no pixel S/N cutoff
+        // XXX probably should be passing the maskVal to this function so we can pass it along here...
+        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) {
             // Can't do anything about it; limp along as best we can
             psErrorClear();
@@ -320 +399 @@
     options->poissonErrorsParams = true;
     options->stats = psStatsAlloc(PSF_STATS);
-    options->radius = maxRadius;
+    options->fitRadius = maxRadius;
+    options->apRadius = maxRadius; // XXX need to decide if aperture mags need a different radius
     options->psfTrendMode = PM_TREND_MAP;
     options->psfTrendNx = xOrder;
@@ -331 +411 @@
 
     pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);
+    pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); // Important for getting a good stack target PSF
 
     pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff);
@@ -343 +424 @@
     pmPSF *psf = psMemIncrRefCounter(try->psf); // Output PSF
     psFree(try);
+
+    if (psTraceGetLevel("psModules.imcombine") >= 1) {
+        psString string = NULL;     // String with values
+        psStringAppend(&string, "Envelope PSF: ");
+        float x = numCols / 2.0, y = numRows / 2.0; // Coordinates of interest
+        for (int j = 4; j < psf->params->n; j++) {
+            pmTrend2D *trend = psf->params->data[j]; // Trend of interest
+            double val = pmTrend2DEval(trend, x, y);
+            double err;
+            switch (trend->mode) {
+              case PM_TREND_POLY_ORD:
+              case PM_TREND_POLY_CHEB:
+                err = NAN;
+                break;
+              case PM_TREND_MAP:
+                err = psImageUnbinPixel(x, y, trend->map->error, trend->map->binning);
+                break;
+              default:
+                psAbort("Unsupported mode: %x", trend->mode);
+            }
+            psStringAppend(&string, "%lf %lf   ", val, err);
+        }
+        psTrace("psModules.imcombine", 1, "%s\n", string);
+        psFree(string);
+    }
 
 #ifdef TESTING
@@ -357 +463 @@
 
         pmReadout *generated = pmReadoutAlloc(NULL); // Generated image
-        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, NULL, NULL, psf, NAN, radius,
+        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, 0, NULL, NULL, psf, NAN, radius,
                                  false, true);
         {