Index: trunk/psphot/src/psPolynomials.c
===================================================================
--- trunk/psphot/src/psPolynomials.c	(revision 4114)
+++ trunk/psphot/src/psPolynomials.c	(revision 4115)
@@ -191,7 +191,8 @@
 // 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,4 +233,5 @@
     // 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,8 +293,9 @@
 // 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,4 +325,5 @@
     // 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,5 +395,5 @@
 }    
 
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
 			    const psVector* x,
 			    const psVector* y,
@@ -413,5 +417,5 @@
     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,5 +428,7 @@
 	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,4 +441,41 @@
 	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);