Changeset 5607 for trunk/psphot/src/pmPSFtry.c

Timestamp:

Nov 25, 2005, 5:00:41 PM (21 years ago)

Author:

eugene

Message:

updating to work with current psModules

File:

• trunk/psphot/src/pmPSFtry.c (modified) (4 diffs)

trunk/psphot/src/pmPSFtry.c

@@ -16 +16 @@
 static void pmPSFtryFree (pmPSFtry *test) {
 
-  if (test == NULL) return;
-
-  psFree (test->psf);
-  psFree (test->sources);
-  psFree (test->modelFLT);
-  psFree (test->modelPSF);
-  psFree (test->metric);
-  psFree (test->fitMag);
-  psFree (test->mask);
-  return;
+    if (test == NULL) return;
+
+    psFree (test->psf);
+    psFree (test->sources);
+    psFree (test->modelFLT);
+    psFree (test->modelPSF);
+    psFree (test->metric);
+    psFree (test->fitMag);
+    psFree (test->mask);
+    return;
 }
 
@@ -95 +95 @@
         // exclude the poor fits
         if (!status) {
-          try->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
-          psFree (model);
-          continue;
+            try->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
+            psFree (model);
+            continue;
         }
         try->modelFLT->data[i] = model;
@@ -136 +136 @@
 
         // otherwise, save the resulting model
-        // pmModelSkyOffset (modelPSF, x, y, RADIUS);
         try->modelPSF->data[i] = modelPSF;
 
@@ -171 +170 @@
 bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
 
-  float dBin;
-  int   nKeep, nSkip;
-
-  // the measured (aperture - fit) magnitudes (dA == try->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
-
-  FILE *f;
-  f = fopen ("apresid.dat", "w");
-  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
-
-  // rflux = ten(0.4*fitMag);
-  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
-  for (int i = 0; i < try->sources->n; i++) {
-    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
-    fprintf (f, "%3d %8.4f %12.5e %8.4f\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i]);
-  }
-  fclose (f);
-
-  // XXX EAM : try 3hi/1lo sigma clipping on the rflux v dap fit
-
-  // find min and max of (1/flux):
-  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-  psVectorStats (stats, rflux, NULL, try->mask, PSFTRY_MASK_ALL);
+    float dBin;
+    int   nKeep, nSkip;
+
+    // the measured (aperture - fit) magnitudes (dA == try->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 (try->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try->sources->n; i++) {
+        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+        rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
+    }
+
+    // find min and max of (1/flux):
+    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+    psVectorStats (stats, rflux, NULL, try->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 (try->sources->n, PS_TYPE_F64);
-    tmp->n = 0;
-
-    // accumulate data within bin range
-    for (int j = 0; j < try->sources->n; j++) {
-      // masked for: bad model fit, outlier in parameters
-      if (try->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
+    // 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 (try->sources->n, PS_TYPE_F64);
+        tmp->n = 0;
+
+        // accumulate data within bin range
+        for (int j = 0; j < try->sources->n; j++) {
+            // masked for: bad model fit, outlier in parameters
+            if (try->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
     
-      // skip points with extreme dA values
-      if (fabs(try->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] = try->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);
-    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;
-
+            // skip points with extreme dA values
+            if (fabs(try->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] = try->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);
+        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 clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
+
+    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);
+
+    try->psf->ApResid = poly->coeff[0];
+    try->psf->dApResid = stats->clippedStdev;
+    try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    psFree (rflux);
+    psFree (rfBin);
+    psFree (daBin);
+    psFree (maskB);
+    psFree (daBinFit);
+    psFree (daResid);
+    psFree (poly);
     psFree (stats);
-    psFree (tmp);
-    psFree (tmp2);
-  }
-
-  // linear clipped fit to rfBin, daBin
-  psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
-  psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-  poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
-
-  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);
-
-  try->psf->ApResid = poly->coeff[0];
-  try->psf->dApResid = stats->clippedStdev;
-  try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
-
-  psFree (rflux);
-  psFree (rfBin);
-  psFree (daBin);
-  psFree (maskB);
-  psFree (daBinFit);
-  psFree (daResid);
-  psFree (poly);
-  psFree (stats);
-  psFree (fitstat);
-
-  return true;
+    psFree (fitstat);
+
+    return true;
 }
 
 bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) {
 
-  // the measured (aperture - fit) magnitudes (dA == try->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 (try->sources->n, PS_TYPE_F64);
-  for (int i = 0; i < try->sources->n; i++) {
-    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
-  }
-
-  // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
-  psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-
-  // XXX EAM 
-  stats->min = 1.0;
-  stats->max = 3.0;
-  stats->clipIter = 3;
-
-  // linear clipped fit to rfBin, daBin
-  psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
-  poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
-  fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
-
-  try->psf->ApResid = poly->coeff[0];
-  try->psf->dApResid = stats->sampleStdev;
-  try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
-
-  FILE *f;
-  f = fopen ("apresid.dat", "w");
-  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
-
-  for (int i = 0; i < try->sources->n; i++) {
-    fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
-  }
-  fclose (f);
-
-  psFree (rflux);
-  psFree (poly);
-  psFree (stats);
-
-  // psFree (daFit);
-  // psFree (daResid);
-
-  return true;
-}
+    // the measured (aperture - fit) magnitudes (dA == try->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 (try->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try->sources->n; i++) {
+        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+        rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
+    }
+
+    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // XXX EAM 
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
+    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    try->psf->ApResid = poly->coeff[0];
+    try->psf->dApResid = stats->sampleStdev;
+    try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    FILE *f;
+    f = fopen ("apresid.dat", "w");
+    if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
+
+    for (int i = 0; i < try->sources->n; i++) {
+        fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
+    }
+    fclose (f);
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}