Changeset 4375 for trunk/psphot/src/pspsf.c

Timestamp:

Jun 24, 2005, 5:38:37 AM (21 years ago)

Author:

eugene

Message:

adding model tests

File:

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

trunk/psphot/src/pspsf.c

@@ -34 +34 @@
     test->modelFLT    = psArrayAlloc (sources->n);
     test->modelPSF    = psArrayAlloc (sources->n);
-    test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
     test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
+    test->ApResid     = 0;
+    test->dApResid    = 0;
+    test->skyBias     = 0;
 
     for (int i = 0; i < test->modelFLT->n; i++) {
+        test->mask->data.U8[i]  = 0;
         test->modelFLT->data[i] = NULL;
         test->modelPSF->data[i] = NULL;
-        test->mask->data.U8[i]  = 0;
+        test->metric->data.F64[i] = 0;
+        test->fitMag->data.F64[i] = 0;
     }   
     return (test);
 }
 
+// test->mask values indicate reason source was rejected:
+// 1: outlier in psf polynomial fit
+// 2: flt model failed to converge
+// 3: psf model failed to converge
+// 4: invalid source photometry
 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
 {
-  
+    bool status;
+    float obsMag;
+    float fitMag;
     float x;
     float y;
-    bool status;
     int Nflt = 0;
     int Npsf = 0;
-    psVector *metricMask = psVectorAlloc (sources->n, PS_TYPE_U8);
+
     pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
 
@@ -66 +77 @@
 
         // set temporary object mask and fit object
-        // fit model as FLT, not PSF (mask & skip poor fits)
+        // fit model as FLT, not PSF
         psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
         status = pmSourceFitModel (source, model, false);
         psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
+        // exclude the poor fits
         if (!status) {
-          test->mask->data.U8[i] = 1;
+          test->mask->data.U8[i] = 2;
           continue;
         }
@@ -84 +96 @@
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
-    
-    // count valid sources
-    int Nkeep = 0;
-    for (int i = 0; i < sources->n; i++) {
-      if (test->mask->data.U8[i]) continue;
-      Nkeep++;
-    }
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
     for (int i = 0; i < test->sources->n; i++) {
+        // masked for: bad model fit, outlier in parameters
+        if (test->mask->data.U8[i]) continue;
+
         psSource *source = test->sources->data[i];
         psModel  *modelFLT = test->modelFLT->data[i];
 
-        // masked for: bad model fit, outlier in parameters
-        if (test->mask->data.U8[i]) continue;
-        if (modelFLT == NULL) {
-          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-          continue;
-        }
-        psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
+        // set shape for this model based on PSF
+        psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); 
         x = source->peak->x;
         y = source->peak->y;
@@ -111 +113 @@
         psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
         status = pmSourceFitModel (source, modelPSF, true);
-        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
         // skip poor fits
         if (!status) {
-          test->mask->data.U8[i] = 1;
-          continue;
+            test->mask->data.U8[i] = 3;
+            goto next_source;
         }
+
+        // otherwise, save the resulting model
         test->modelPSF->data[i] = modelPSF;
+
+        // XXX : use a different aperture radius from the fit radius?
+        // XXX : use a different estimator for the local sky?
+        // XXX : pass 'source' as input?
+        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+            test->mask->data.U8[i] = 4;
+            goto next_source;
+        }           
+
+        test->metric->data.F64[i] = obsMag - fitMag;
+        test->fitMag->data.F64[i] = fitMag;
         Npsf ++;
+
+    next_source:
+        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
     }
     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
@@ -125 +143 @@
     DumpModelFits (test->modelPSF, "modelsPSF.dat");
 
-    // stage 4: measure PSF metric (Ap-Fit)
-    // should I be saving the aperture photometry for all sources?
-    // save the metric for all stars
-    for (int i = 0; i < test->sources->n; i++) {
-
-        float obsSum = 0;
-        float fitSum = 0;
-      
-        // is this metricMask redundant with test->mask ?
-        metricMask->data.U8[i] = 1;
-        test->metric->data.F64[i] = 0;
-
-        if (test->mask->data.U8[i]) continue;
-
-        psModel *model = test->modelPSF->data[i];
-        if (model == NULL) {
-          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-          continue;
-        }
-
-        psImage *image = ((psSource *)test->sources->data[i])->pixels;
-        psImage *mask  = ((psSource *)test->sources->data[i])->mask;
-
-        // this metric is Ap-Fit
-        float sky = model->params->data.F32[0];
-        for (int ix = 0; ix < image->numCols; ix++) {
-            for (int iy = 0; iy < image->numRows; iy++) {
-                if (mask->data.U8[iy][ix]) continue;
-                obsSum += image->data.F32[iy][ix] - sky;
-                fitSum += psModelEval (model, image, ix, iy) - sky;
-            }
-        }
-        // fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
-        // keep the good metrics
-        if ((fitSum > 0) && (obsSum > 0)) {
-            metricMask->data.U8[i] = 0;
-            test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
-        } else {
-          test->mask->data.U8[i] = 1;
-        }
-    }
-
-    // XXX use robust stats, not clipped stats
-    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    psVectorStats (stats, test->metric, NULL, metricMask, 1);
-    test->metricStats = stats;
-    psLogMsg ("psphot.pspsf", 3, "test model %s, metric: %f +/- %f\n", modelName, stats->clippedMean, stats->clippedStdev);
+    // XXX this function wants aperture radius from pmSourcePhotometry
+    pmPSFMetricModel (test, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
+              modelName, test->ApResid, test->dApResid, test->skyBias);
     return (test);
 }
@@ -181 +156 @@
     // construct the fit vectors from the collection of objects
     // use the mask to ignore missing fits 
-    psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
 
@@ -208 +183 @@
             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
         }
-        // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
 
         psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
+        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
         // psPolynomial2DDump (psf->params->data[i]);
-
-        // count valid sources
-        int Nkeep = 0;
-        for (int j = 0; j < mask->n; j++) {
-          if (mask->data.U8[j]) continue;
-          Nkeep++;
-        }
-        psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
-
     }
     return (true);
@@ -241 +207 @@
 }
 
+bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+
+  float dBin;
+  int   nKeep, nSkip;
+  
+
+  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+  //     dA = dAo + dsky/flux
+  //   where flux is the flux of the star
+  // we fit this trend to find the infinite flux aperture correction (dAo),
+  //   the nominal sky bias (dsky), and the error on dAo
+  // the values of dA are contaminated by stars with close neighbors in the aperture
+  //   we use an outlier rejection to avoid this bias
+
+  // rflux = ten(0.4*fitMag);
+  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < test->sources->n; i++) {
+    if (test->mask->data.U8[i]) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  }
+  // psScalar *t1 = psScalar(0.4);
+  // psVector *v1 = psBinaryOp (NULL, test->fitMag, "*", t1);
+  // psVector *rflux = psUnaryOp (NULL, v1, "ten");
+  // psFree (t1);
+  // psFree (v1);
+
+  // find min and max of (1/flux):
+  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
+  
+  // build binned versions of rflux, metric
+  dBin = (stats->max - stats->min) / 10.0;
+  psVector *rfBin = psVectorCreate (stats->min, stats->max, dBin, PS_TYPE_F64);
+  psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
+  psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
+  psFree (stats);
+
+  psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
+
+  for (int i = 0; i < daBin->n; i++) {
+
+    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    tmp->n = 0;
+
+    // accumulate data within bin range
+    for (int j = 0; j < test->sources->n; j++) {
+      // masked for: bad model fit, outlier in parameters
+      if (test->mask->data.U8[j]) continue;
+    
+      // skip points with extreme dA values
+      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+
+      // skip points outside of this bin
+      if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
+      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
+
+      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->n ++;
+    }
+
+    // is this a valid point?
+    maskB->data.U8[i] = 0;
+    if (tmp->n < 2) {
+      maskB->data.U8[i] = 1;
+      continue;
+    } 
+
+    // dA values are contaminated with low outliers 
+    // measure statistics only on upper 50% of points
+    psVectorSort (tmp, tmp);
+    nKeep = 0.5*tmp->n;
+    nSkip = tmp->n - nKeep;
+
+    psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
+    for (int j = 0; j < tmp2->n; j++) {
+      tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
+    }
+
+    stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, tmp2, NULL, NULL, 0);
+    psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+    daBin->data.F64[i] = stats->sampleMedian;
+  }
+
+  // linear fit to rfBin, daBin
+  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
+
+  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
+
+  stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+  test->ApResid = poly->coeff[0];
+  test->dApResid = stats->clippedStdev;
+  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+  return true;
+}