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

Timestamp:

Jun 5, 2005, 8:41:28 PM (21 years ago)

Author:

eugene

Message:

using source masks to select PSF stars

File:

• trunk/psphot/src/psPolynomials.c (modified) (1 diff)

trunk/psphot/src/psPolynomials.c

@@ -453 +453 @@
                             const psVector* dz)
 {
+    PS_VECTOR_CHECK_NULL(mask, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_NULL(dz, NULL);
+
     psVector *zFit   = NULL;
-    psVector *zResid = NULL;
-    psStats  *stats  = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+    psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
         poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
         zFit   = Polynomial2DEvalVectorD (poly, x, y);
-        zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
-
-        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
         stats  = psVectorStats (stats, zResid, NULL, mask, 1);
-        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", 
-                 stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
+        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+                 stats->sampleMean, stats->sampleStdev);
 
         // set mask if pts are not valid
-        int n = 0;
+        // 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 (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
-              mask->data.U8[i] = 1;
-              continue;
-            } else {
-              mask->data.U8[i] = 0;
-            }         
-            n++;
-        }
-    }
+          if (mask->data.U8[i]) continue;
+          if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+            mask->data.U8[i] = 1;
+            continue;
+          }       
+        }
+        psFree (zFit);
+    }
+    psFree (zResid);
+    psFree (stats);
     return (poly);
 }