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

Timestamp:

Jun 4, 2005, 5:22:13 PM (21 years ago)

Author:

eugene

Message:

running sources in brightness order, using a mask

File:

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

trunk/psphot/src/psPolynomials.c

@@ -191 +191 @@
 // XXX EAM : test version of 1d fitting
 psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* yErr)
+                                             psVector *mask,
+                                             const psVector *x,
+                                             const psVector *y,
+                                             const psVector *yErr)
 {
     // I think this is 1 dimension down
@@ -232 +233 @@
     // Build the B and A data structs.
     for (int k = 0; k < x->n; k++) {
+        if ((mask != NULL) && mask->U8[k]) continue;
         xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
       
@@ -291 +293 @@
 // XXX EAM : test version of 2d fitting
 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* z,
-        const psVector* zErr)
+                                             psVector* mask,
+                                             const psVector* x,
+                                             const psVector* y,
+                                             const psVector* z,
+                                             const psVector* zErr)
 {
     // I think this is 1 dimension down
@@ -322 +325 @@
     // Build the B and A data structs.
     for (int k  = 0; k < x->n; k++) {
+      if ((mask != NULL) && mask->data.U8[i]) continue;
         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
       
@@ -391 +395 @@
 }    
 
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
                             const psVector* x,
                             const psVector* y,
@@ -413 +417 @@
     for (int N = 0; N < 3; N++) {
         // XXX EAM : this would be better defined with an element mask
-        poly   = VectorFitPolynomial2DOrd_EAM (poly, X, Y, Z, dZ);
+        poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
         zFit   = Polynomial2DEvalVectorD (poly, x, y);
         zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
@@ -424 +428 @@
         int n = 0;
         for (int i = 0; i < zResid->n; i++) {
-            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) continue;
+            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
+              continue;
+            } 
             X->data.F64[n]  =  x->data.F64[i];
             Y->data.F64[n]  =  y->data.F64[i];
@@ -435 +441 @@
         Z->n = n;
         dZ->n = n;
+    }
+    return (poly);
+}
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+                            psVector* mask,
+                            const psVector* x,
+                            const psVector* y,
+                            const psVector* z,
+                            const psVector* dz)
+{
+    psVector *zFit   = NULL;
+    psVector *zResid = NULL;
+    psStats  *stats  = NULL;
+
+    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);
+        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);
+
+        // set mask if pts are not valid
+        int n = 0;
+        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++;
+        }
     }
     return (poly);