Index: trunk/psphot/src/psPolynomials.c
===================================================================
--- trunk/psphot/src/psPolynomials.c	(revision 4115)
+++ trunk/psphot/src/psPolynomials.c	(revision 4116)
@@ -453,30 +453,37 @@
 			    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);
 }