Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 9540)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 9625)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-13 22:04:58 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-10-18 18:37:38 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -42,4 +42,5 @@
 
 # define USE_GAUSS_JORDAN 1
+# define USE_ROBUST_STATS_FOR_CLIPPING 1
 
 #define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
@@ -811,4 +812,5 @@
     // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
     stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -859,8 +861,24 @@
 
         stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
-        psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
+        # if (USE_ROBUST_STATS_FOR_CLIPPING)
+            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->robustMedian;
+        # else
+
+            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->sampleMedian;
+        # endif
 
         // set mask if pts are not valid
@@ -872,6 +890,5 @@
             }
 
-            if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
-                    (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
+            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -1246,4 +1263,5 @@
     // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
     stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -1252,8 +1270,8 @@
         psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n", N);
         psS32 Nkeep = 0;
-        if (psTraceGetLevel("psLib.math") >= 6) {
+        if (psTraceGetLevel("psLib.math") >= 7) {
             if (mask != NULL) {
                 for (psS32 i = 0 ; i < mask->n ; i++) {
-                    psTrace("psLib.math", 6,  "mask[%d] is %d\n", i, mask->data.U8[i]);
+                    psTrace("psLib.math", 7,  "mask[%d] is %d\n", i, mask->data.U8[i]);
                 }
             }
@@ -1282,9 +1300,9 @@
         }
 
-        if (psTraceGetLevel("psLib.math") >= 6) {
+        if (psTraceGetLevel("psLib.math") >= 7) {
             if (mask != NULL) {
                 for (psS32 i = 0 ; i < mask->n ; i++) {
                     if (!((mask != NULL) && (mask->data.U8[i] & maskValue))) {
-                        psTrace("psLib.math", 6,  "(f, fit)[%d] is (%f, %f).  resid is (%f)\n",
+                        psTrace("psLib.math", 7,  "(f, fit)[%d] is (%f, %f).  resid is (%f)\n",
                                 i, f->data.F32[i], fit->data.F32[i], resid->data.F64[i]);
                     }
@@ -1300,8 +1318,18 @@
             return(NULL);
         }
-        psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
+        # if (USE_ROBUST_STATS_FOR_CLIPPING)
+            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->robustMedian;
+        # else
+
+            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->sampleMedian;
+        # endif
 
         // set mask if pts are not valid
@@ -1313,6 +1341,5 @@
             }
 
-            if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
-                    (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
+            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
                 if (fit->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -1331,5 +1358,5 @@
         }
 
-        psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        psTrace("psLib.math", 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
         psFree(fit);
     }
@@ -1766,4 +1793,5 @@
     // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
     stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -1820,8 +1848,18 @@
         }
 
-        psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
+        # if (USE_ROBUST_STATS_FOR_CLIPPING)
+            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->robustMedian;
+        # else
+
+            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->sampleMedian;
+        # endif
 
         // set mask if pts are not valid
@@ -1833,6 +1871,5 @@
             }
 
-            if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
-                    (resid->data.F64[i] - stats->sampleMedian < minClipValue))  {
+            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue))  {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
@@ -2334,4 +2371,5 @@
     // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
     stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
     psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
@@ -2388,8 +2426,18 @@
             return(NULL);
         }
-        psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
-        psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
-        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
-        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
+        # if (USE_ROBUST_STATS_FOR_CLIPPING)
+            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->robustMedian;
+        # else
+
+            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
+        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
+        psF32 minClipValue = -minClipSigma*stats->robustStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
+        psF32 clipMedian = stats->sampleMedian;
+        # endif
 
         // set mask if pts are not valid
@@ -2401,6 +2449,5 @@
             }
 
-            if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
-                    (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
+            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
                 if (f->type.type == PS_TYPE_F64) {
                     psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",