Changeset 5958 for branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c

Timestamp:

Jan 9, 2006, 6:46:03 PM (21 years ago)

Author:

magnier

Message:

updated objects code with ApTrend concept (was in eam_rel9_b0)

File:

• branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c (modified) (6 diffs)

branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -119 +119 @@
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
 
-    // make this optional?
-    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
-
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
@@ -169 +166 @@
 
     }
+    psfTry->psf->nPSFstars = Npsf;
+
     psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
@@ -177 +176 @@
     // XXX this function wants aperture radius for pmSourcePhotometry
     pmPSFtryMetric (psfTry, RADIUS);
-    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
-              modelName,
-              psfTry->psf->ApResid,
-              psfTry->psf->dApResid,
-              psfTry->psf->skyBias);
-
+    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
+              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
     return (psfTry);
 }
 
-
 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
 {
-
-    float dBin;
-    int   nKeep, nSkip;
-
     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
@@ -210 +200 @@
     }
 
-    // find min and max of (1/flux):
-    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-    psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
-
-    // build binned versions of rflux, metric
-    dBin = (stats->max - stats->min) / 10.0;
-    psVector *rfBin = psVectorCreate(NULL, 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);
-
-    // group data in daBin bins, measure lower 50% mean
-    for (int i = 0; i < daBin->n; i++) {
-
-        psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
-        tmp->n = 0;
-
-        // accumulate data within bin range
-        for (int j = 0; j < psfTry->sources->n; j++) {
-            // masked for: bad model fit, outlier in parameters
-            if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
-                continue;
-
-            // skip points with extreme dA values
-            if (fabs(psfTry->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] = psfTry->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;
-            psFree (tmp);
-            continue;
-        }
-
-        // dA values are contaminated with low outliers
-        // measure statistics only on upper 50% of points
-        // this would be easier if we could sort in reverse:
-        //
-        // psVectorSort (tmp, tmp);
-        // tmp->n = 0.5*tmp->n;
-        // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        // psVectorStats (stats, tmp, NULL, NULL, 0);
-        // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
-
-        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;
-
-        psFree (stats);
-        psFree (tmp);
-        psFree (tmp2);
-    }
-
-    // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
-
-    // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
-
-    // XXX EAM : replace this when the above version is implemented
-    // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
-
-    psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
-    psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
-
-    stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
-    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
-
-    psfTry->psf->ApResid = poly->coeff[0];
-    psfTry->psf->dApResid = stats->clippedStdev;
-    psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+    // use 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // fit ApResid only to rflux, ignore x,y variations for now
+    // linear clipped fit of ApResid to rflux
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, rflux);
+    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    psfTry->psf->ApTrend->coeff[0][0][0] = poly->coeff[0];
+    psfTry->psf->ApTrend->coeff[0][0][1] = 0;
+
+    psfTry->psf->skyBias  = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+    psfTry->psf->ApResid  = poly->coeff[0];
+    psfTry->psf->dApResid = stats->sampleStdev;
 
     psFree (rflux);
-    psFree (rfBin);
-    psFree (daBin);
-    psFree (maskB);
-    psFree (daBinFit);
-    psFree (daResid);
     psFree (poly);
     psFree (stats);
@@ -317 +225 @@
     return true;
 }
+
+/*
+  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+*/
+