Changeset 5068 for trunk/psphot/src/psPolynomials.c

Timestamp:

Sep 19, 2005, 12:26:40 PM (21 years ago)

Author:

eugene

Message:

fixing errors related to apResid and skyBias

File:

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

trunk/psphot/src/psPolynomials.c

@@ -497 +497 @@
 
 // XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial1D* VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly,
+                                              psStats *stats,
+                                              psVector* mask,
+                                              psMaskType maskValue,
+                                              const psVector* z,
+                                              const psVector* dz,
+                                              const psVector* x)
+{
+    psVector *zFit   = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < 3; N++) {
+        int Nkeep = 0;
+
+        poly   = VectorFitPolynomial1D_EAM (poly, mask, maskValue, z, dz, x);
+        zFit   = Polynomial1DEvalVector_EAM (poly, x);
+        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
+        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
+        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+                 stats->sampleMedian, stats->sampleStdev);
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range 
+        // recovery is not allowed with this scheme
+        for (int i = 0; i < zResid->n; i++) {
+            if (mask->data.U8[i]) continue;
+            if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 2*stats->sampleStdev) {
+                mask->data.U8[i] |= 0x01;
+                continue;
+            }       
+            Nkeep ++;
+        }
+
+        psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", 
+                 Nkeep, x->n);
+
+        psFree (zFit);
+    }
+    psFree (zResid);
+    return (poly);
+}
+
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
 psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
                                               psStats *stats,
@@ -510 +558 @@
 
     // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
-    stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
+        int Nkeep = 0;
+
         poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
         zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
@@ -519 +569 @@
         stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
-                 stats->sampleMean, stats->sampleStdev);
+                 stats->sampleMedian, stats->sampleStdev);
 
         // set mask if pts are not valid
@@ -526 +576 @@
         for (int i = 0; i < zResid->n; i++) {
             if (mask->data.U8[i]) continue;
-            if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
-                mask->data.U8[i] &= 0x01;
+            if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 3*stats->sampleStdev) {
+                mask->data.U8[i] |= 0x01;
                 continue;
             }       
-        }
+            Nkeep ++;
+        }
+
+        psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", 
+                 Nkeep, x->n);
+
+//      psTrace (".psphot.RobustFit", 4, "model pars: %f %f %f %f\n", 
+//               poly->coeff[0][0], poly->coeff[1][0], 
+//               poly->coeff[0][1], poly->coeff[1][1]);
+
+        psTrace (".psphot.RobustFit", 4, "model pars: %f\n", poly->coeff[0][0]);
+
         psFree (zFit);
     }