Index: trunk/psLib/src/math/psPolynomialMD.c
===================================================================
--- trunk/psLib/src/math/psPolynomialMD.c	(revision 17075)
+++ trunk/psLib/src/math/psPolynomialMD.c	(revision 17147)
@@ -153,4 +153,34 @@
 }
 
+// Calculate the standard deviation of the fit
+static void polynomialMDStdev(psPolynomialMD *poly, // Polynomial for which to measure stdev
+                              psVector *deviations, // Deviations, or NULL
+                              const psArray *coords, // Array of coordinates
+                              const psVector *values, // Measured values
+                              const psVector *mask // Mask for values
+                              )
+{
+    assert(poly);
+    assert(coords);
+    assert(values);
+
+    double rms = 0.0;                   // Root mean square deviation
+    int numValues = values->n;          // Number of values
+    int numGood = numValues;            // Number of good values
+    for (int i = 0; i < numValues; i++) {
+        if (mask && mask->data.U8[i]) {
+            numGood--;
+            continue;
+        }
+        float diff = values->data.F32[i] - psPolynomialMDEval(poly, coords->data[i]);
+        if (deviations) {
+            deviations->data.F32[i] = diff;
+        }
+        rms += PS_SQR(diff);
+    }
+    poly->stdevFit = sqrt(rms / (double)numGood);
+
+    return;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -195,4 +225,6 @@
     poly->dim = dim;
     poly->orders = psVectorCopy(NULL, orders, PS_TYPE_U8);
+    poly->numFit = 0;
+    poly->stdevFit = NAN;
 
     int numTerms = 1;                   // Number of terms in polynomial
@@ -301,4 +333,6 @@
         return false;
     }
+
+    polynomialMDStdev(poly, NULL, coordsArray, values, NULL);
 
     return true;
@@ -389,19 +423,9 @@
         }
 
-        // Evaluate distribution
-        double rms = 0.0;               // Root mean square deviation
-        for (int i = 0; i < numValues; i++) {
-            if (ownMask->data.U8[i]) {
-                continue;
-            }
-            float diff = values->data.F32[i] - psPolynomialMDEval(poly, coordsArray->data[i]);
-            deviations->data.F32[i] = diff;
-            rms += PS_SQR(diff);
-        }
-        rms = sqrt(rms / (double)numGood);
-        psTrace("psLib.math", 7, "RMS from %d points is %lf\n", numGood, rms);
+        polynomialMDStdev(poly, deviations, coordsArray, values, ownMask);
+        psTrace("psLib.math", 7, "RMS from %d points is %lf\n", numGood, poly->stdevFit);
 
         // Reject
-        float limit = rej * rms;        // Rejection limit
+        float limit = rej * poly->stdevFit; // Rejection limit
         for (int i = 0; i < numValues; i++) {
             if (ownMask->data.U8[i]) {
@@ -433,4 +457,5 @@
             return false;
         }
+        polynomialMDStdev(poly, NULL, coordsArray, values, ownMask);
     }
 
@@ -446,2 +471,3 @@
     return true;
 }
+