Index: /trunk/psLib/src/dataManip/psFunctions.c
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.c	(revision 888)
+++ /trunk/psLib/src/dataManip/psFunctions.c	(revision 889)
@@ -6,6 +6,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-04 01:03:11 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -71,4 +71,6 @@
  routines for generating both uniformly distributed numbers and the
  Gaussian distribution as well.
+ 
+ GUS: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -82,4 +84,5 @@
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
+    gauss->n = Npts;
     gsl_rng_env_setup();
     T = gsl_rng_default;
@@ -120,7 +123,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (float *) psAlloc(nY * sizeof(float));
@@ -129,4 +132,5 @@
     }
 
+
     return(newPoly);
 }
@@ -143,11 +147,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float));
-        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -173,15 +177,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float));
-            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -410,7 +414,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (double *) psAlloc(nY * sizeof(double));
@@ -433,11 +437,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double));
-        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (double *) psAlloc(nZ * sizeof(double));
@@ -463,15 +467,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double));
-            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (double *) psAlloc(nZ * sizeof(double));
Index: /trunk/psLib/src/dataManip/psFunctions.h
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.h	(revision 888)
+++ /trunk/psLib/src/dataManip/psFunctions.h	(revision 889)
@@ -18,4 +18,9 @@
           );
 
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts);
+
 /*****************************************************************************/
 
@@ -205,5 +210,5 @@
 /** Evaluate 1D polynomial (double precision) */
 double
-psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+psDEvalPolynomial1D(double x,  ///< Value at which to evaluate
                     const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
                    );
@@ -211,5 +216,5 @@
 /** Evaluate 2D polynomial (double precision) */
 double
-psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial2D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
@@ -218,5 +223,5 @@
 /** Evaluate 3D polynomial (double precision) */
 double
-psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial3D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     double z,  ///< Value z at which to evaluate
@@ -226,5 +231,5 @@
 /** Evaluate 4D polynomial (double precision) */
 double
-psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+psDEvalPolynomial4D(double w,  ///< Value w at which to evaluate
                     double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 888)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 889)
@@ -85,6 +85,12 @@
     float binSize = 0.0;
 
+    if (n == 0) {
+        psAbort(__func__, "psHistogram requested with 0 bins");
+    }
+
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32);
+    newHist->bounds->n = newHist->bounds->nalloc;
+
     binSize = (upper - lower) / (float) n;
     for (i=0;i<n+1;i++) {
@@ -92,4 +98,8 @@
     }
     newHist->nums =  psVectorAlloc(n, PS_TYPE_S32);
+    newHist->nums->n = newHist->nums->nalloc;
+    for (i=0;i<newHist->nums->n;i++) {
+        newHist->nums->data.S32[i] = 0;
+    }
     newHist->minNum = 0;
     newHist->maxNum = 0;
@@ -99,5 +109,6 @@
 }
 
-
+// When this is called, the number of data elements in bounds
+// is n.  Therefore, the number of bins is n-1.
 psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
 {
@@ -107,8 +118,13 @@
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
+    newHist->bounds->n = newHist->bounds->nalloc;
     for (i=0;i<bounds->n;i++) {
         newHist->bounds->data.F32[i] = bounds->data.F32[i];
     }
     newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32);
+    newHist->nums->n = newHist->nums->nalloc;
+    for (i=0;i<newHist->nums->n;i++) {
+        newHist->nums->data.S32[i] = 0;
+    }
 
     newHist->minNum = 0;
@@ -119,6 +135,8 @@
 }
 
-void psHistogramFree(psHistogram *restrict myHist)
-{
+void psHistogramFree(psHistogram *myHist)
+{
+    psVectorFree(myHist->bounds);
+    psVectorFree(myHist->nums);
     psFree(myHist);
 }
@@ -150,15 +168,51 @@
 
     numBins = out->nums->n;
-    for (i=0;i<in->n;i++) {
-        // Check if this pixel is masked, and if so, skip it.
-        if (!(mask->data.S32[i] & maskVal)) {
-            // Check if this pixel is below the minimum value, and if so
-            // count it, then skip it.
+
+    if (mask != NULL) {
+        for (i=0;i<in->n;i++) {
+            // Check if this pixel is masked, and if so, skip it.
+            if (!(mask->data.U8[i] & maskVal)) {
+                // Check if this pixel is below the minimum value, and if so
+                // count it, then skip it.
+                if (in->data.F32[i] < out->bounds->data.F32[0]) {
+                    out->minNum++;
+
+                    // Check if this pixel is above the maximum value, and if so
+                    // count it, then skip it.
+                } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
+                    out->maxNum++;
+                } else {
+                    // If this is a uniform histogram, determining the correct
+                    // number is trivial.
+                    if (out->uniform == true) {
+                        binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
+
+                        binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
+                                        binSize);
+                        (out->nums->data.S32[binNum])++;
+                        // If this is a non-uniform histogram, determining the correct
+                        // bin number requires a bit more work.
+                    } else {
+                        // GUS: This is slow.  Put a smarter algorithm here to
+                        // find the correct bin number (bin search, probably)
+                        for (j=0;j<(out->bounds->n)-1;j++) {
+                            if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
+                                    (in->data.F32[i] <= out->bounds->data.S32[j+1])) {
+                                (out->nums->data.S32[j])++;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    } else {
+        for (i=0;i<in->n;i++) {
             if (in->data.F32[i] < out->bounds->data.F32[0]) {
+                // Check if this pixel is below the minimum value, and if so
+                // count it, then skip it.
                 out->minNum++;
-
+            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
                 // Check if this pixel is above the maximum value, and if so
                 // count it, then skip it.
-            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
                 out->maxNum++;
             } else {
@@ -171,4 +225,5 @@
                                     binSize);
                     (out->nums->data.S32[binNum])++;
+
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
@@ -185,4 +240,5 @@
             }
         }
+
     }
     return(out);
@@ -471,4 +527,5 @@
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     unsortedVector->n = unsortedVector->nalloc;
+
     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     sortedVector->n = sortedVector->nalloc;
@@ -494,5 +551,5 @@
                 if ((rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -504,10 +561,10 @@
             for (i=0;i<myVector->n;i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                unsortedVector->data.F32[i] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->data.F32[i] = myVector->data.F32[i];
             }
         }
@@ -596,25 +653,27 @@
     float rangeMin = 0.0;
     float rangeMax = 0.0;
-    psStats *stats2 = NULL;
 
     // Determine if the number of data points exceed a threshold which will
     // cause to generate robust stats, as opposed to exact stats.
-    if (myVector->n > stats->sampleLimit) {
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleUQ = stats2->robustUQ;
-        stats->sampleLQ = stats2->robustLQ;
-
-        // Free temporary data buffers.
-        psStatsFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
+    /* GUS: When IfA provides the algorithm, insert is here.
+        if (myVector->n > stats->sampleLimit) {
+            psStats *stats2 = NULL;
+            // Calculate the robust quartiles.
+            stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
+            p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
+     
+            // Store the robust quartiles into the sample quartile members.
+            stats->sampleUQ = stats2->robustUQ;
+            stats->sampleLQ = stats2->robustLQ;
+     
+            // Free temporary data buffers.
+            psStatsFree(stats2);
+     
+            // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
+            stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
+     
+            return;
+        }
+    */
 
     // Determine how many data points fit inside this min/max range
@@ -640,5 +699,5 @@
                         (rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -647,5 +706,5 @@
                 if ((rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -657,10 +716,10 @@
             for (i=0;i<myVector->n;i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                unsortedVector->data.F32[i] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->data.F32[i] = myVector->data.F32[i];
             }
         }
@@ -679,4 +738,5 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
+    // GUS: This is the
 }
 
@@ -691,5 +751,5 @@
         PS_STAT_ROBUST_QUARTILE
     I have included all that computation in a single function, as opposed to
-    breaking it across several functions for one primary reason: the all
+    breaking it across several functions for one primary reason: they all
     require the same basic initial processing steps (calculate the histogram,
     etc.) however there is no currently defined means, in the SDRS, to
@@ -719,8 +779,7 @@
     // is really required.
 
-    if (0.0 == stats->sampleUQ) {
-
-        // GUS: fix this
-        //        stats->options = stats->options | PS_STAT_SAMPLE_UQ;
+    if (isnan(stats->sampleUQ) ||
+            isnan(stats->sampleLQ)) {
+        stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
         p_psVectorSampleQuartiles(myVector,
                                   maskVector,
@@ -729,13 +788,4 @@
     }
 
-    if (isnan(stats->sampleLQ)) {
-        // GUS: fix this
-        //        stats->options = stats->options | PS_STAT_SAMPLE_LQ;
-        p_psVectorSampleQuartiles(myVector,
-                                  maskVector,
-                                  maskVal,
-                                  stats);
-    }
-
     // Compute the initial bin size of the robust histogram.
     sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f;
@@ -753,6 +803,8 @@
     robustHistogram = psHistogramAlloc(stats->min,
                                        stats->max,
-                                       binSize);
-    // Populate the histogram arrat.
+                                       1 + (int) ((stats->max - stats->min) / binSize));
+
+
+    // Populate the histogram array.
     // GUS: fix this
     //    robustHistogram = psHistogramVector(robustHistogram, myVector);
@@ -811,4 +863,7 @@
         stats->robustLQ = 0.0;
     }
+    stats->robustNfit = 0.0;
+    stats->robustN50 = 0.0;
+    psHistogramFree(robustHistogram);
 }
 
@@ -836,5 +891,5 @@
 
     if (0 != isnan(stats->sampleMean)) {
-        psAbort(__func__, "stats->sampleMean == NAN");
+        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     }
     mean = stats->sampleMean;
@@ -920,12 +975,17 @@
     }
 
-    tmpMask = psVectorAlloc(myVector->nalloc, maskVector->type.type);
-
-    tmpMask->n = maskVector->n;
-    for (i=0;i<tmpMask->n;i++) {
-        tmpMask->data.U8[i] = maskVector->data.U8[i];
+
+    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+    tmpMask->n = myVector->n;
+
+    if (maskVector != NULL) {
+        for (i=0;i<tmpMask->n;i++) {
+            tmpMask->data.U8[i] = maskVector->data.U8[i];
+        }
     }
 
     // 1. Compute the sample median.
+    // GUS: This seems odd.  Verify with IfA that we want to calculate the
+    // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
 
@@ -942,7 +1002,9 @@
 
     for (i=0;i<stats->clipIter;i++) {
+        //        printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i,
+        //                clippedMean, clippedStdev);
         for (j=0;j<myVector->n;j++) {
             // a) Exclude all values x_i for which |x_i - x| > K * stdev
-            if ( fabs(myVector->data.F32[j] - clippedMean) >
+            if (fabs(myVector->data.F32[j] - clippedMean) >
                     (stats->clipSigma * clippedStdev)) {
                 tmpMask->data.U8[i] = 0xff;
@@ -994,6 +1056,8 @@
     if (in->type.type != PS_TYPE_F32) {
         psAbort(__func__,
-                "Only data type PS_TYPE_F32 is currently supported.");
-    }
+                "Only data type PS_TYPE_F32 is currently supported (0x%x).",
+                in->type.type);
+    }
+
     if (mask != NULL) {
         if (in->n != mask->n) {
@@ -1018,5 +1082,8 @@
 
     // ************************************************************************
+    // GUS: The Stdev calculation requires the mean.  Should we assume the
+    // mean has already been calculated?  Or should we always calculate it?
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
+        p_psVectorSampleMean(in, mask, maskVal, stats);
         p_psVectorSampleStdev(in, mask, maskVal, stats);
     }
Index: /trunk/psLib/src/dataManip/psStats.h
===================================================================
--- /trunk/psLib/src/dataManip/psStats.h	(revision 888)
+++ /trunk/psLib/src/dataManip/psStats.h	(revision 889)
@@ -81,5 +81,5 @@
 typedef struct
 {
-    const psVector *restrict bounds; ///< Bounds for the bins (type F32)
+    psVector *bounds;                ///< Bounds for the bins (type F32)
     psVector *nums;                  ///< Number in each of the bins (INT)
     int minNum;                      ///< Number below the minimum
@@ -102,5 +102,5 @@
 
 /** Destructor \ingroup MathGroup **/
-void psHistogramFree(psHistogram *restrict myHist);          ///< Histogram to destroy
+void psHistogramFree(psHistogram *myHist);          ///< Histogram to destroy
 
 
Index: /trunk/psLib/src/image/psImage.c
===================================================================
--- /trunk/psLib/src/image/psImage.c	(revision 888)
+++ /trunk/psLib/src/image/psImage.c	(revision 889)
@@ -9,8 +9,10 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-04 03:15:47 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  That is the routine used to generate matrices.
  */
 
Index: /trunk/psLib/src/image/psImageStats.c
===================================================================
--- /trunk/psLib/src/image/psImageStats.c	(revision 888)
+++ /trunk/psLib/src/image/psImageStats.c	(revision 889)
@@ -19,4 +19,9 @@
 
 /// This routine must determine the various statistics for the image.
+/*****************************************************************************
+    GUS: verify that image/mask have the
+ correct types
+ sizes
+ *****************************************************************************/
 psStats *psImageStats(psStats *stats,
                       psImage *in,
@@ -26,15 +31,39 @@
     psVector *junkData=NULL;
     psVector *junkMask=NULL;
+    int ptr = 0;
+    int i = 0;
+    int j = 0;
 
     junkData = psVectorAlloc(in->numRows * in->numCols, in->type.type);
-    junkMask = psVectorAlloc(mask->numRows * mask->numCols, mask->type.type);
-
-    // GUS: figure out mask types
-    junkData->data.F32 = (float *) in->data.F32;
-    junkMask->data.F32 = (float *) mask->data.F32;
-    stats = psVectorStats(stats, junkData, junkMask, maskVal);
+    junkData->n = junkData->nalloc;
+    ptr=0;
+    for (i=0;i<in->numRows;i++) {
+        for (j=0;j<in->numCols;j++) {
+            junkData->data.F32[ptr++] = in->data.F32[i][j];
+        }
+    }
+
+    if (mask != NULL) {
+        // GUS: verify that mask data is PS_TYPE_U8.
+        // GUS: figure out mask types
+
+        junkMask = psVectorAlloc(mask->numRows * mask->numCols,
+                                 mask->type.type);
+        junkMask->n = junkMask->nalloc;
+        // GUS: Is there a more efficient way to do this?
+        ptr=0;
+        for (i=0;i<mask->numRows;i++) {
+            for (j=0;j<mask->numCols;j++) {
+                junkMask->data.U8[ptr++] = mask->data.U8[i][j];
+            }
+        }
+
+        stats = psVectorStats(stats, junkData, junkMask, maskVal);
+        psVectorFree(junkMask);
+    } else {
+        stats = psVectorStats(stats, junkData, NULL, 0);
+    }
 
     psVectorFree(junkData);
-    psVectorFree(junkMask);
 
     return(stats);
@@ -44,4 +73,8 @@
     NOTE: We assume that the psHistogram structure out has already been
     allocated and initialized.
+ 
+    GUS: verify that image/mask have the
+ correct types
+ sizes
  *****************************************************************************/
 psHistogram *psImageHistogram(psHistogram *out,
@@ -52,16 +85,45 @@
     psVector *junkData=NULL;
     psVector *junkMask=NULL;
+    int ptr = 0;
+    int i = 0;
+    int j = 0;
 
     junkData = psVectorAlloc(in->numRows * in->numCols, in->type.type);
-    junkMask = psVectorAlloc(mask->numRows * mask->numCols, mask->type.type);
-
-    // GUS: figure out mask types
-    junkData->data.F32 = (float *) in->data.F32;
-    junkMask->data.F32 = (float *) mask->data.F32;
-
-    out = psHistogramVector(out, junkData, junkMask, maskVal);
-
+    junkData->n = junkData->nalloc;
+
+    // GUS: Is there a more efficient way to do this?  memcopy() won't work.
+    // Can we trick the junkData structure to use the image buffer, then
+    // untrick it before we deallocate it (so we won't deallocate that buffer
+    // twice?
+
+    // GUS: Make sure you have the numRows/NumCols in the right place.
+    ptr=0;
+    for (i=0;i<in->numRows;i++) {
+        for (j=0;j<in->numCols;j++) {
+            junkData->data.F32[ptr++] = in->data.F32[i][j];
+        }
+    }
+
+    if (mask != NULL) {
+        // GUS: verify that mask data is PS_TYPE_U8.
+        // GUS: figure out mask types
+
+        junkMask = psVectorAlloc(mask->numRows * mask->numCols,
+                                 mask->type.type);
+        junkMask->n = junkMask->nalloc;
+
+        // GUS: Is there a more efficient way to do this?
+        ptr=0;
+        for (i=0;i<mask->numRows;i++) {
+            for (j=0;j<mask->numCols;j++) {
+                junkMask->data.U8[ptr++] = mask->data.U8[i][j];
+            }
+        }
+        out = psHistogramVector(out, junkData, junkMask, maskVal);
+        psVectorFree(junkMask);
+    } else {
+        out = psHistogramVector(out, junkData, NULL, 0);
+    }
     psVectorFree(junkData);
-    psVectorFree(junkMask);
 
     return(out);
@@ -91,6 +153,8 @@
     int maxChebyPoly = 0;
 
+
     // Create the sums[][] data structure.  This will hold the LHS of equation
     // 29 in the ADD: sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
+
     sums = (float **) psAlloc(coeffs->nX * sizeof(float *));
     for (i=0;i<coeffs->nX;i++) {
@@ -105,5 +169,5 @@
     // Determine how many Chebyshev polynomials are needed, then create them.
     maxChebyPoly = coeffs->nX;
-    if (coeffs->nX > coeffs->nY) {
+    if (coeffs->nY > coeffs->nX) {
         maxChebyPoly = coeffs->nY;
     }
@@ -133,4 +197,5 @@
                                  psEvalPolynomial1D((float) x, chebPolys[i]) *
                                  psEvalPolynomial1D((float) y, chebPolys[j]);
+
                 }
             }
@@ -160,9 +225,15 @@
     }
 
-
     // Free the Chebyshev polynomials that were created in this routine.
     for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
+        psPolynomial1DFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    // Free some data
+    for (i=0;i<coeffs->nX;i++) {
+        psFree(sums[i]);
+    }
+    psFree(sums);
 
     return(coeffs);
@@ -199,5 +270,5 @@
     // Determine how many Chebyshev polynomials are needed, then create them.
     maxChebyPoly = coeffs->nX;
-    if (coeffs->nX > coeffs->nY) {
+    if (coeffs->nY > coeffs->nX) {
         maxChebyPoly = coeffs->nY;
     }
@@ -219,5 +290,4 @@
         }
     }
-
 
     for (x=0;x<input->numRows;x++) {
@@ -227,7 +297,6 @@
                 for (j=0;j<coeffs->nY;j++) {
                     polySum+= psEvalPolynomial1D((float)x, chebPolys[i]) *
-                              psEvalPolynomial1D((float)y, chebPolys[y]) *
+                              psEvalPolynomial1D((float)y, chebPolys[j]) *
                               coeffs->coeff[i][j];
-
                 }
             }
@@ -236,9 +305,15 @@
     }
 
-
     // Free the Chebyshev polynomials that were created in this routine.
     for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
+        psPolynomial1DFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    // Free some data
+    for (i=0;i<coeffs->nX;i++) {
+        psFree(sums[i]);
+    }
+    psFree(sums);
 
     return(0);
Index: /trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- /trunk/psLib/src/imageops/psImageStats.c	(revision 888)
+++ /trunk/psLib/src/imageops/psImageStats.c	(revision 889)
@@ -19,4 +19,9 @@
 
 /// This routine must determine the various statistics for the image.
+/*****************************************************************************
+    GUS: verify that image/mask have the
+ correct types
+ sizes
+ *****************************************************************************/
 psStats *psImageStats(psStats *stats,
                       psImage *in,
@@ -26,15 +31,39 @@
     psVector *junkData=NULL;
     psVector *junkMask=NULL;
+    int ptr = 0;
+    int i = 0;
+    int j = 0;
 
     junkData = psVectorAlloc(in->numRows * in->numCols, in->type.type);
-    junkMask = psVectorAlloc(mask->numRows * mask->numCols, mask->type.type);
-
-    // GUS: figure out mask types
-    junkData->data.F32 = (float *) in->data.F32;
-    junkMask->data.F32 = (float *) mask->data.F32;
-    stats = psVectorStats(stats, junkData, junkMask, maskVal);
+    junkData->n = junkData->nalloc;
+    ptr=0;
+    for (i=0;i<in->numRows;i++) {
+        for (j=0;j<in->numCols;j++) {
+            junkData->data.F32[ptr++] = in->data.F32[i][j];
+        }
+    }
+
+    if (mask != NULL) {
+        // GUS: verify that mask data is PS_TYPE_U8.
+        // GUS: figure out mask types
+
+        junkMask = psVectorAlloc(mask->numRows * mask->numCols,
+                                 mask->type.type);
+        junkMask->n = junkMask->nalloc;
+        // GUS: Is there a more efficient way to do this?
+        ptr=0;
+        for (i=0;i<mask->numRows;i++) {
+            for (j=0;j<mask->numCols;j++) {
+                junkMask->data.U8[ptr++] = mask->data.U8[i][j];
+            }
+        }
+
+        stats = psVectorStats(stats, junkData, junkMask, maskVal);
+        psVectorFree(junkMask);
+    } else {
+        stats = psVectorStats(stats, junkData, NULL, 0);
+    }
 
     psVectorFree(junkData);
-    psVectorFree(junkMask);
 
     return(stats);
@@ -44,4 +73,8 @@
     NOTE: We assume that the psHistogram structure out has already been
     allocated and initialized.
+ 
+    GUS: verify that image/mask have the
+ correct types
+ sizes
  *****************************************************************************/
 psHistogram *psImageHistogram(psHistogram *out,
@@ -52,16 +85,45 @@
     psVector *junkData=NULL;
     psVector *junkMask=NULL;
+    int ptr = 0;
+    int i = 0;
+    int j = 0;
 
     junkData = psVectorAlloc(in->numRows * in->numCols, in->type.type);
-    junkMask = psVectorAlloc(mask->numRows * mask->numCols, mask->type.type);
-
-    // GUS: figure out mask types
-    junkData->data.F32 = (float *) in->data.F32;
-    junkMask->data.F32 = (float *) mask->data.F32;
-
-    out = psHistogramVector(out, junkData, junkMask, maskVal);
-
+    junkData->n = junkData->nalloc;
+
+    // GUS: Is there a more efficient way to do this?  memcopy() won't work.
+    // Can we trick the junkData structure to use the image buffer, then
+    // untrick it before we deallocate it (so we won't deallocate that buffer
+    // twice?
+
+    // GUS: Make sure you have the numRows/NumCols in the right place.
+    ptr=0;
+    for (i=0;i<in->numRows;i++) {
+        for (j=0;j<in->numCols;j++) {
+            junkData->data.F32[ptr++] = in->data.F32[i][j];
+        }
+    }
+
+    if (mask != NULL) {
+        // GUS: verify that mask data is PS_TYPE_U8.
+        // GUS: figure out mask types
+
+        junkMask = psVectorAlloc(mask->numRows * mask->numCols,
+                                 mask->type.type);
+        junkMask->n = junkMask->nalloc;
+
+        // GUS: Is there a more efficient way to do this?
+        ptr=0;
+        for (i=0;i<mask->numRows;i++) {
+            for (j=0;j<mask->numCols;j++) {
+                junkMask->data.U8[ptr++] = mask->data.U8[i][j];
+            }
+        }
+        out = psHistogramVector(out, junkData, junkMask, maskVal);
+        psVectorFree(junkMask);
+    } else {
+        out = psHistogramVector(out, junkData, NULL, 0);
+    }
     psVectorFree(junkData);
-    psVectorFree(junkMask);
 
     return(out);
@@ -91,6 +153,8 @@
     int maxChebyPoly = 0;
 
+
     // Create the sums[][] data structure.  This will hold the LHS of equation
     // 29 in the ADD: sums[k][l] = SUM { image(x,y) * Tk(x) * Tl(y) }
+
     sums = (float **) psAlloc(coeffs->nX * sizeof(float *));
     for (i=0;i<coeffs->nX;i++) {
@@ -105,5 +169,5 @@
     // Determine how many Chebyshev polynomials are needed, then create them.
     maxChebyPoly = coeffs->nX;
-    if (coeffs->nX > coeffs->nY) {
+    if (coeffs->nY > coeffs->nX) {
         maxChebyPoly = coeffs->nY;
     }
@@ -133,4 +197,5 @@
                                  psEvalPolynomial1D((float) x, chebPolys[i]) *
                                  psEvalPolynomial1D((float) y, chebPolys[j]);
+
                 }
             }
@@ -160,9 +225,15 @@
     }
 
-
     // Free the Chebyshev polynomials that were created in this routine.
     for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
+        psPolynomial1DFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    // Free some data
+    for (i=0;i<coeffs->nX;i++) {
+        psFree(sums[i]);
+    }
+    psFree(sums);
 
     return(coeffs);
@@ -199,5 +270,5 @@
     // Determine how many Chebyshev polynomials are needed, then create them.
     maxChebyPoly = coeffs->nX;
-    if (coeffs->nX > coeffs->nY) {
+    if (coeffs->nY > coeffs->nX) {
         maxChebyPoly = coeffs->nY;
     }
@@ -219,5 +290,4 @@
         }
     }
-
 
     for (x=0;x<input->numRows;x++) {
@@ -227,7 +297,6 @@
                 for (j=0;j<coeffs->nY;j++) {
                     polySum+= psEvalPolynomial1D((float)x, chebPolys[i]) *
-                              psEvalPolynomial1D((float)y, chebPolys[y]) *
+                              psEvalPolynomial1D((float)y, chebPolys[j]) *
                               coeffs->coeff[i][j];
-
                 }
             }
@@ -236,9 +305,15 @@
     }
 
-
     // Free the Chebyshev polynomials that were created in this routine.
     for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
+        psPolynomial1DFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    // Free some data
+    for (i=0;i<coeffs->nX;i++) {
+        psFree(sums[i]);
+    }
+    psFree(sums);
 
     return(0);
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 888)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 889)
@@ -6,6 +6,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-04 01:03:11 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -71,4 +71,6 @@
  routines for generating both uniformly distributed numbers and the
  Gaussian distribution as well.
+ 
+ GUS: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -82,4 +84,5 @@
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
+    gauss->n = Npts;
     gsl_rng_env_setup();
     T = gsl_rng_default;
@@ -120,7 +123,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (float *) psAlloc(nY * sizeof(float));
@@ -129,4 +132,5 @@
     }
 
+
     return(newPoly);
 }
@@ -143,11 +147,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float));
-        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -173,15 +177,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float));
-            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -410,7 +414,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (double *) psAlloc(nY * sizeof(double));
@@ -433,11 +437,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double));
-        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (double *) psAlloc(nZ * sizeof(double));
@@ -463,15 +467,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double));
-            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (double *) psAlloc(nZ * sizeof(double));
Index: /trunk/psLib/src/math/psPolynomial.h
===================================================================
--- /trunk/psLib/src/math/psPolynomial.h	(revision 888)
+++ /trunk/psLib/src/math/psPolynomial.h	(revision 889)
@@ -18,4 +18,9 @@
           );
 
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts);
+
 /*****************************************************************************/
 
@@ -205,5 +210,5 @@
 /** Evaluate 1D polynomial (double precision) */
 double
-psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+psDEvalPolynomial1D(double x,  ///< Value at which to evaluate
                     const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
                    );
@@ -211,5 +216,5 @@
 /** Evaluate 2D polynomial (double precision) */
 double
-psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial2D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
@@ -218,5 +223,5 @@
 /** Evaluate 3D polynomial (double precision) */
 double
-psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial3D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     double z,  ///< Value z at which to evaluate
@@ -226,5 +231,5 @@
 /** Evaluate 4D polynomial (double precision) */
 double
-psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+psDEvalPolynomial4D(double w,  ///< Value w at which to evaluate
                     double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 888)
+++ /trunk/psLib/src/math/psSpline.c	(revision 889)
@@ -6,6 +6,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-04 01:03:11 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -71,4 +71,6 @@
  routines for generating both uniformly distributed numbers and the
  Gaussian distribution as well.
+ 
+ GUS: There is no way to seed the random generator.
  *****************************************************************************/
 psVector *psGaussianDev(float mean,
@@ -82,4 +84,5 @@
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
+    gauss->n = Npts;
     gsl_rng_env_setup();
     T = gsl_rng_default;
@@ -120,7 +123,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (float *) psAlloc(nY * sizeof(float));
@@ -129,4 +132,5 @@
     }
 
+
     return(newPoly);
 }
@@ -143,11 +147,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float));
-        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -173,15 +177,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float));
-            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (float *) psAlloc(nZ * sizeof(float));
@@ -410,7 +414,7 @@
     newPoly->nY = nY;
 
-    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double *));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char *));
     for (x=0;x<nX;x++) {
         newPoly->coeff[x]    = (double *) psAlloc(nY * sizeof(double));
@@ -433,11 +437,11 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
-    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double **));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char **));
     for (x=0;x<nX;x++) {
-        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double));
-        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double));
-        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char));
+        newPoly->coeff[x]    = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double *));
+        newPoly->mask[x]     = (char **)  psAlloc(nY * sizeof(char *));
         for (y=0;y<nY;y++) {
             newPoly->coeff[x][y]    = (double *) psAlloc(nZ * sizeof(double));
@@ -463,15 +467,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
-    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double ***));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char ***));
     for (w=0;w<nW;w++) {
-        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double));
-        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char));
+        newPoly->coeff[w]    = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double **));
+        newPoly->mask[w]     = (char ***)  psAlloc(nX * sizeof(char **));
         for (x=0;x<nX;x++) {
-            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double));
-            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double));
-            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char));
+            newPoly->coeff[w][x]    = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double *));
+            newPoly->mask[w][x]     = (char **) psAlloc(nY * sizeof(char *));
             for (y=0;y<nY;y++) {
                 newPoly->coeff[w][x][y]    = (double *) psAlloc(nZ * sizeof(double));
Index: /trunk/psLib/src/math/psSpline.h
===================================================================
--- /trunk/psLib/src/math/psSpline.h	(revision 888)
+++ /trunk/psLib/src/math/psSpline.h	(revision 889)
@@ -18,4 +18,9 @@
           );
 
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts);
+
 /*****************************************************************************/
 
@@ -205,5 +210,5 @@
 /** Evaluate 1D polynomial (double precision) */
 double
-psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+psDEvalPolynomial1D(double x,  ///< Value at which to evaluate
                     const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
                    );
@@ -211,5 +216,5 @@
 /** Evaluate 2D polynomial (double precision) */
 double
-psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial2D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
@@ -218,5 +223,5 @@
 /** Evaluate 3D polynomial (double precision) */
 double
-psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+psDEvalPolynomial3D(double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
                     double z,  ///< Value z at which to evaluate
@@ -226,5 +231,5 @@
 /** Evaluate 4D polynomial (double precision) */
 double
-psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+psDEvalPolynomial4D(double w,  ///< Value w at which to evaluate
                     double x,  ///< Value x at which to evaluate
                     double y,  ///< Value y at which to evaluate
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 888)
+++ /trunk/psLib/src/math/psStats.c	(revision 889)
@@ -85,6 +85,12 @@
     float binSize = 0.0;
 
+    if (n == 0) {
+        psAbort(__func__, "psHistogram requested with 0 bins");
+    }
+
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32);
+    newHist->bounds->n = newHist->bounds->nalloc;
+
     binSize = (upper - lower) / (float) n;
     for (i=0;i<n+1;i++) {
@@ -92,4 +98,8 @@
     }
     newHist->nums =  psVectorAlloc(n, PS_TYPE_S32);
+    newHist->nums->n = newHist->nums->nalloc;
+    for (i=0;i<newHist->nums->n;i++) {
+        newHist->nums->data.S32[i] = 0;
+    }
     newHist->minNum = 0;
     newHist->maxNum = 0;
@@ -99,5 +109,6 @@
 }
 
-
+// When this is called, the number of data elements in bounds
+// is n.  Therefore, the number of bins is n-1.
 psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
 {
@@ -107,8 +118,13 @@
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
+    newHist->bounds->n = newHist->bounds->nalloc;
     for (i=0;i<bounds->n;i++) {
         newHist->bounds->data.F32[i] = bounds->data.F32[i];
     }
     newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32);
+    newHist->nums->n = newHist->nums->nalloc;
+    for (i=0;i<newHist->nums->n;i++) {
+        newHist->nums->data.S32[i] = 0;
+    }
 
     newHist->minNum = 0;
@@ -119,6 +135,8 @@
 }
 
-void psHistogramFree(psHistogram *restrict myHist)
-{
+void psHistogramFree(psHistogram *myHist)
+{
+    psVectorFree(myHist->bounds);
+    psVectorFree(myHist->nums);
     psFree(myHist);
 }
@@ -150,15 +168,51 @@
 
     numBins = out->nums->n;
-    for (i=0;i<in->n;i++) {
-        // Check if this pixel is masked, and if so, skip it.
-        if (!(mask->data.S32[i] & maskVal)) {
-            // Check if this pixel is below the minimum value, and if so
-            // count it, then skip it.
+
+    if (mask != NULL) {
+        for (i=0;i<in->n;i++) {
+            // Check if this pixel is masked, and if so, skip it.
+            if (!(mask->data.U8[i] & maskVal)) {
+                // Check if this pixel is below the minimum value, and if so
+                // count it, then skip it.
+                if (in->data.F32[i] < out->bounds->data.F32[0]) {
+                    out->minNum++;
+
+                    // Check if this pixel is above the maximum value, and if so
+                    // count it, then skip it.
+                } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
+                    out->maxNum++;
+                } else {
+                    // If this is a uniform histogram, determining the correct
+                    // number is trivial.
+                    if (out->uniform == true) {
+                        binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
+
+                        binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
+                                        binSize);
+                        (out->nums->data.S32[binNum])++;
+                        // If this is a non-uniform histogram, determining the correct
+                        // bin number requires a bit more work.
+                    } else {
+                        // GUS: This is slow.  Put a smarter algorithm here to
+                        // find the correct bin number (bin search, probably)
+                        for (j=0;j<(out->bounds->n)-1;j++) {
+                            if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
+                                    (in->data.F32[i] <= out->bounds->data.S32[j+1])) {
+                                (out->nums->data.S32[j])++;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    } else {
+        for (i=0;i<in->n;i++) {
             if (in->data.F32[i] < out->bounds->data.F32[0]) {
+                // Check if this pixel is below the minimum value, and if so
+                // count it, then skip it.
                 out->minNum++;
-
+            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
                 // Check if this pixel is above the maximum value, and if so
                 // count it, then skip it.
-            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
                 out->maxNum++;
             } else {
@@ -171,4 +225,5 @@
                                     binSize);
                     (out->nums->data.S32[binNum])++;
+
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
@@ -185,4 +240,5 @@
             }
         }
+
     }
     return(out);
@@ -471,4 +527,5 @@
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     unsortedVector->n = unsortedVector->nalloc;
+
     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     sortedVector->n = sortedVector->nalloc;
@@ -494,5 +551,5 @@
                 if ((rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -504,10 +561,10 @@
             for (i=0;i<myVector->n;i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                unsortedVector->data.F32[i] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->data.F32[i] = myVector->data.F32[i];
             }
         }
@@ -596,25 +653,27 @@
     float rangeMin = 0.0;
     float rangeMax = 0.0;
-    psStats *stats2 = NULL;
 
     // Determine if the number of data points exceed a threshold which will
     // cause to generate robust stats, as opposed to exact stats.
-    if (myVector->n > stats->sampleLimit) {
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleUQ = stats2->robustUQ;
-        stats->sampleLQ = stats2->robustLQ;
-
-        // Free temporary data buffers.
-        psStatsFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
+    /* GUS: When IfA provides the algorithm, insert is here.
+        if (myVector->n > stats->sampleLimit) {
+            psStats *stats2 = NULL;
+            // Calculate the robust quartiles.
+            stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
+            p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
+     
+            // Store the robust quartiles into the sample quartile members.
+            stats->sampleUQ = stats2->robustUQ;
+            stats->sampleLQ = stats2->robustLQ;
+     
+            // Free temporary data buffers.
+            psStatsFree(stats2);
+     
+            // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
+            stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
+     
+            return;
+        }
+    */
 
     // Determine how many data points fit inside this min/max range
@@ -640,5 +699,5 @@
                         (rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -647,5 +706,5 @@
                 if ((rangeMin <= myVector->data.F32[i]) &&
                         (myVector->data.F32[i] <= rangeMax)) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
             }
@@ -657,10 +716,10 @@
             for (i=0;i<myVector->n;i++) {
                 if (!(maskVal & maskVector->data.U8[i])) {
-                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                unsortedVector->data.F32[i] = maskVector->data.F32[i];
+                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->data.F32[i] = myVector->data.F32[i];
             }
         }
@@ -679,4 +738,5 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
+    // GUS: This is the
 }
 
@@ -691,5 +751,5 @@
         PS_STAT_ROBUST_QUARTILE
     I have included all that computation in a single function, as opposed to
-    breaking it across several functions for one primary reason: the all
+    breaking it across several functions for one primary reason: they all
     require the same basic initial processing steps (calculate the histogram,
     etc.) however there is no currently defined means, in the SDRS, to
@@ -719,8 +779,7 @@
     // is really required.
 
-    if (0.0 == stats->sampleUQ) {
-
-        // GUS: fix this
-        //        stats->options = stats->options | PS_STAT_SAMPLE_UQ;
+    if (isnan(stats->sampleUQ) ||
+            isnan(stats->sampleLQ)) {
+        stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
         p_psVectorSampleQuartiles(myVector,
                                   maskVector,
@@ -729,13 +788,4 @@
     }
 
-    if (isnan(stats->sampleLQ)) {
-        // GUS: fix this
-        //        stats->options = stats->options | PS_STAT_SAMPLE_LQ;
-        p_psVectorSampleQuartiles(myVector,
-                                  maskVector,
-                                  maskVal,
-                                  stats);
-    }
-
     // Compute the initial bin size of the robust histogram.
     sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f;
@@ -753,6 +803,8 @@
     robustHistogram = psHistogramAlloc(stats->min,
                                        stats->max,
-                                       binSize);
-    // Populate the histogram arrat.
+                                       1 + (int) ((stats->max - stats->min) / binSize));
+
+
+    // Populate the histogram array.
     // GUS: fix this
     //    robustHistogram = psHistogramVector(robustHistogram, myVector);
@@ -811,4 +863,7 @@
         stats->robustLQ = 0.0;
     }
+    stats->robustNfit = 0.0;
+    stats->robustN50 = 0.0;
+    psHistogramFree(robustHistogram);
 }
 
@@ -836,5 +891,5 @@
 
     if (0 != isnan(stats->sampleMean)) {
-        psAbort(__func__, "stats->sampleMean == NAN");
+        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     }
     mean = stats->sampleMean;
@@ -920,12 +975,17 @@
     }
 
-    tmpMask = psVectorAlloc(myVector->nalloc, maskVector->type.type);
-
-    tmpMask->n = maskVector->n;
-    for (i=0;i<tmpMask->n;i++) {
-        tmpMask->data.U8[i] = maskVector->data.U8[i];
+
+    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+    tmpMask->n = myVector->n;
+
+    if (maskVector != NULL) {
+        for (i=0;i<tmpMask->n;i++) {
+            tmpMask->data.U8[i] = maskVector->data.U8[i];
+        }
     }
 
     // 1. Compute the sample median.
+    // GUS: This seems odd.  Verify with IfA that we want to calculate the
+    // median here, not the mean.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
 
@@ -942,7 +1002,9 @@
 
     for (i=0;i<stats->clipIter;i++) {
+        //        printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i,
+        //                clippedMean, clippedStdev);
         for (j=0;j<myVector->n;j++) {
             // a) Exclude all values x_i for which |x_i - x| > K * stdev
-            if ( fabs(myVector->data.F32[j] - clippedMean) >
+            if (fabs(myVector->data.F32[j] - clippedMean) >
                     (stats->clipSigma * clippedStdev)) {
                 tmpMask->data.U8[i] = 0xff;
@@ -994,6 +1056,8 @@
     if (in->type.type != PS_TYPE_F32) {
         psAbort(__func__,
-                "Only data type PS_TYPE_F32 is currently supported.");
-    }
+                "Only data type PS_TYPE_F32 is currently supported (0x%x).",
+                in->type.type);
+    }
+
     if (mask != NULL) {
         if (in->n != mask->n) {
@@ -1018,5 +1082,8 @@
 
     // ************************************************************************
+    // GUS: The Stdev calculation requires the mean.  Should we assume the
+    // mean has already been calculated?  Or should we always calculate it?
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
+        p_psVectorSampleMean(in, mask, maskVal, stats);
         p_psVectorSampleStdev(in, mask, maskVal, stats);
     }
Index: /trunk/psLib/src/math/psStats.h
===================================================================
--- /trunk/psLib/src/math/psStats.h	(revision 888)
+++ /trunk/psLib/src/math/psStats.h	(revision 889)
@@ -81,5 +81,5 @@
 typedef struct
 {
-    const psVector *restrict bounds; ///< Bounds for the bins (type F32)
+    psVector *bounds;                ///< Bounds for the bins (type F32)
     psVector *nums;                  ///< Number in each of the bins (INT)
     int minNum;                      ///< Number below the minimum
@@ -102,5 +102,5 @@
 
 /** Destructor \ingroup MathGroup **/
-void psHistogramFree(psHistogram *restrict myHist);          ///< Histogram to destroy
+void psHistogramFree(psHistogram *myHist);          ///< Histogram to destroy
 
 
Index: /trunk/psLib/src/mathtypes/psImage.c
===================================================================
--- /trunk/psLib/src/mathtypes/psImage.c	(revision 888)
+++ /trunk/psLib/src/mathtypes/psImage.c	(revision 889)
@@ -9,8 +9,10 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-04 03:15:47 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  That is the routine used to generate matrices.
  */
 
Index: /trunk/psLib/src/sys/psLogMsg.c
===================================================================
--- /trunk/psLib/src/sys/psLogMsg.c	(revision 888)
+++ /trunk/psLib/src/sys/psLogMsg.c	(revision 889)
@@ -29,4 +29,5 @@
 static int logDest = PS_LOG_TO_STDERR; // where to log messages
 static int logLevel = PS_LOG_INFO;     // log all messages at this or above
+//static FILE *p_logFP = NULL;
 // The following variables control
 // which information is displayed in
Index: /trunk/psLib/src/sysUtils/psHash.c
===================================================================
--- /trunk/psLib/src/sysUtils/psHash.c	(revision 888)
+++ /trunk/psLib/src/sysUtils/psHash.c	(revision 889)
@@ -1,7 +1,18 @@
-/******************************************************************************
-    An interface to hash tables for Pan-STARRS
- 
-    Collisions are handled by simple linked lists
- *****************************************************************************/
+/** @file  psHash.c
+ *
+ *  @brief Contains support for basic hashing functions.
+ *
+ *  This file will hold the functions for defining a hash table with arbitrary
+ *  data types, allocating/deallocating that has table, adding and removing
+ *  data from that hash table, and listing all keys defined in the hash table.
+ *
+ *  @author Robert Lupton, Princeton University
+ *  @author George Gusciora, MHPCC
+ *   
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
 #include <stdlib.h>
 #include <stdio.h>
@@ -12,4 +23,5 @@
 #include "psString.h"
 #include "psTrace.h"
+#include "psAbort.h"
 
 psDlist *psHashKeyList(psHash *table)
@@ -142,11 +154,27 @@
     // bucket.
     psHashBucket *optr = NULL;
+    char *tmpchar = NULL;
 
     // NOTE: this is NOT a good hash function!  See Knuth.
     //
-    for (int i = 0, len = strlen(key); i < len; i++) {
-        hash = (hash << 1) ^ key[i];
-    }
-    hash &= (table->nbucket - 1);
+    //    for (int i = 0, len = strlen(key); i < len; i++) {
+    //        hash = (hash << 1) ^ key[i];
+    //    }
+    //    hash &= (table->nbucket - 1);
+
+    // This hash algorithm is from Sedgewick.  NOTE: must reread to ensure that
+    // the size of the hash table is not required to be a prime number.
+    tmpchar = (char *) key;
+    for (hash = 0; *tmpchar != '\0'; tmpchar++) {
+        hash = (64 * hash + *tmpchar) % (table->nbucket);
+    }
+
+    // NOTE: This should not be necessary, but for now, I'm checking bounds
+    // anyway.
+    if ((hash < 0) ||
+            (hash >= table->nbucket)) {
+        psAbort(__func__,"Internal hash function out of range (%d)",
+                hash);
+    }
 
     ptr = table->buckets[hash];
Index: /trunk/psLib/src/sysUtils/psHash.h
===================================================================
--- /trunk/psLib/src/sysUtils/psHash.h	(revision 888)
+++ /trunk/psLib/src/sysUtils/psHash.h	(revision 889)
@@ -1,2 +1,19 @@
+/** @file  psHash.c
+ *
+ *  @brief Contains support for basic hashing functions.
+ *
+ *  This file will hold the prototypes for defining a hash table with arbitrary
+ *  data types, allocating/deallocating that has table, adding and removing
+ *  data from that hash table, and listing all keys defined in the hash table.
+ *
+ *  @author Robert Lupton, Princeton University
+ *  @author George Gusciora, MHPCC
+ *   
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-07 00:32:53 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
 #if !defined(PS_HASH_H)
 #define PS_HASH_H
Index: /trunk/psLib/src/sysUtils/psLogMsg.c
===================================================================
--- /trunk/psLib/src/sysUtils/psLogMsg.c	(revision 888)
+++ /trunk/psLib/src/sysUtils/psLogMsg.c	(revision 889)
@@ -29,4 +29,5 @@
 static int logDest = PS_LOG_TO_STDERR; // where to log messages
 static int logLevel = PS_LOG_INFO;     // log all messages at this or above
+//static FILE *p_logFP = NULL;
 // The following variables control
 // which information is displayed in