Index: trunk/psLib/src/math/psClip.c
===================================================================
--- trunk/psLib/src/math/psClip.c	(revision 11619)
+++ trunk/psLib/src/math/psClip.c	(revision 11755)
@@ -15,5 +15,5 @@
 #include "psClip.h"
 
-// No-op; purpose of function is to identify the type
+// No-op; purpose of function is merely to identify the type
 static void clipParamsFree(psClipParams *params)
 {
@@ -44,5 +44,77 @@
 
 
-long psClip(psClipParams *params, psVector *values, psVector *mask, const psVector *errors)
+long psClipMinMax(const psClipParams *params, const psVector *values, psVector *mask)
+{
+    PS_ASSERT_VECTOR_NON_NULL(values, -1);
+    PS_ASSERT_VECTOR_NON_NULL(mask, -1);
+    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_MASK, -1);
+    PS_ASSERT_PTR(params, -1);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, -1);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracLow, 0.0, -1);
+    PS_ASSERT_FLOAT_LESS_THAN(params->fracLow, 1.0, -1);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracHigh, 0.0, -1);
+    PS_ASSERT_FLOAT_LESS_THAN(params->fracHigh, 1.0, -1);
+    PS_ASSERT_INT_NONZERO(params->clipped, -1);
+
+    if (params->fracLow == 0.0 && params->fracHigh == 0.0) {
+        // No min-max rejection desired
+        return 0;
+    }
+
+    psMaskType masked = params->masked; // Indicates masked values
+    psMaskType clipped = params->clipped; // Indicates clipped values
+    masked |= clipped;                  // Make sure we're also masking clipped values
+    psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask
+    long totalMasked = 0;               // Total number of pixels masked
+
+    // Apply fracLow,fracHigh if there are enough values
+
+    // Run through to get number of operational values
+    long numValid = 0;                  // Number of valid values
+    for (long i = 0; i < mask->n; i++) {
+        if (!(maskData[i] & masked)) {
+            numValid++;
+        }
+    }
+    psTrace("psLib.math", 2, "%ld valid values.\n", numValid);
+
+    // XXX: Not sure if sorting provides the fastest implementation.  It might be quicker to do a linear pass
+    // through the data, pulling out the highest M and lowest N values,
+
+    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep
+    if (numValid * keepFrac >= params->numKeep) {
+        psTrace("psLib.math", 1, "Applying min/max clipping.\n");
+        psVector *index = psVectorSortIndex(NULL, values);
+        int numLow = numValid * params->fracLow; // Number of low values to clip
+        int numHigh = numValid * params->fracHigh; // Number of high values to clip
+        // Low values
+        psS32 *indexData = index->data.S32; // Dereference index
+        long numMasked = 0;             // Number masked
+        for (long i = 0; i < index->n && numMasked < numLow; i++) {
+            // Don't count the ones that are already masked
+            if (! (maskData[indexData[i]] & masked)) {
+                maskData[indexData[i]] |= clipped;
+                numMasked++;
+            }
+        }
+        totalMasked += numMasked;
+        numMasked = 0;
+        // High values
+        for (long i = values->n - 1;  i >= 0 && numMasked < numHigh; i--) {
+            // Don't count the ones that are already masked
+            if (! (maskData[indexData[i]] & masked)) {
+                maskData[indexData[i]] |= clipped;
+                numMasked++;
+            }
+        }
+        totalMasked += numMasked;
+        psFree(index);
+    }
+
+    return totalMasked;
+}
+
+
+long psClipReject(psClipParams *params, const psVector *values, psVector *mask, const psVector *errors)
 {
     PS_ASSERT_VECTOR_NON_NULL(values, -1);
@@ -57,5 +129,7 @@
         PS_ASSERT_VECTOR_TYPE_EQUAL(values, errors, -1);
     }
-    if (params->iter > 0 && (!isfinite(params->rej) || params->rej <= 0)) {
+    PS_ASSERT_INT_NONZERO(params->meanStat, -1);
+    PS_ASSERT_INT_NONZERO(params->stdevStat, -1);
+    if (!isfinite(params->rej) || params->rej <= 0) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Rejection limit (%f) is not valid.\n", params->rej);
         return -1;
@@ -65,54 +139,6 @@
     psMaskType clipped = params->clipped; // Indicates clipped values
     masked |= clipped;                  // Make sure we're also masking clipped values
-    psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask
-    long totalMasked = 0;               // Total number of pixels masked
 
-    // Immediate min/max rejection: apply fracLow,fracHigh if there are enough values
-    if (params->fracLow != 0.0 || params->fracHigh != 0.0) {
-        // Run through to get number of operational values
-        long numValid = 0;                  // Number of valid values
-        for (long i = 0; i < mask->n; i++) {
-            if (!(maskData[i] & masked)) {
-                numValid++;
-            }
-        }
-        psTrace("psLib.math", 2, "%ld valid values.\n", numValid);
-
-        float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep
-        if (numValid * keepFrac >= params->numKeep) {
-            psTrace("psLib.math", 1, "Applying min/max clipping.\n");
-            psVector *index = psVectorSortIndex(NULL, values);
-            int numLow = numValid * params->fracLow; // Number of low values to clip
-            int numHigh = numValid * params->fracHigh; // Number of high values to clip
-            // Low values
-            psS32 *indexData = index->data.S32; // Dereference index
-            long numMasked = 0;             // Number masked
-            for (long i = 0; i < index->n && numMasked < numLow; i++) {
-                // Don't count the ones that are already masked
-                if (! (maskData[indexData[i]] & masked)) {
-                    maskData[indexData[i]] |= clipped;
-                    numMasked++;
-                }
-            }
-            totalMasked += numMasked;
-            numMasked = 0;
-            // High values
-            for (long i = values->n - 1;  i >= 0 && numMasked < numHigh; i--) {
-                // Don't count the ones that are already masked
-                if (! (maskData[indexData[i]] & masked)) {
-                    maskData[indexData[i]] |= clipped;
-                    numMasked++;
-                }
-            }
-            totalMasked += numMasked;
-            psFree(index);
-
-            // Turn off min/max rejection so that future calls will not trigger it
-            params->fracHigh = 0.0;
-            params->fracLow = 0.0;
-        }
-    }
-
-    // Iterate, with rejection
+    // Get statistics
     psStats *stats = psStatsAlloc(params->meanStat | params->stdevStat);
     if (!psVectorStats(stats, values, errors, mask, masked)) {
@@ -124,38 +150,29 @@
     params->stdev = psStatsGetValue(stats, params->stdevStat);
 
-
     #define REJECT_CASE(TYPE) \
 case PS_TYPE_##TYPE: { \
         ps##TYPE *valuesData = values->data.TYPE; /* Dereference for speed */ \
-        for (int i = 0; i < params->iter; i++) { \
-            if (errors) { \
-                ps##TYPE *errorsData = errors->data.TYPE; \
-                for (int j = 0; j < values->n; j++) { \
-                    if (!(maskData[j] & masked) && \
-                            fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \
-                        maskData[j] |= clipped; \
-                        totalMasked++; \
-                    } \
-                } \
-            } else { \
-                ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \
-                for (int j = 0; j < values->n; j++) { \
-                    if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \
-                        maskData[j] |= clipped; \
-                        totalMasked++; \
-                    } \
+        psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; /* Dereference mask for speed */ \
+        if (errors) { \
+            ps##TYPE *errorsData = errors->data.TYPE; \
+            for (int j = 0; j < values->n; j++) { \
+                if (!(maskData[j] & masked) && \
+                        fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \
+                    maskData[j] |= clipped; \
+                    totalMasked++; \
                 } \
             } \
-            \
-            if (!psVectorStats(stats, values, errors, mask, masked)) { \
-                psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on values.\n"); \
-                psFree(stats); \
-                return -1; \
+        } else { \
+            ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \
+            for (int j = 0; j < values->n; j++) { \
+                if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \
+                    maskData[j] |= clipped; \
+                    totalMasked++; \
+                } \
             } \
-            params->mean = psStatsGetValue(stats, params->meanStat); \
-            params->stdev = psStatsGetValue(stats, params->stdevStat); \
         } \
     }
 
+    long totalMasked = 0;               // Total number of pixels masked
     switch (values->type.type) {
         REJECT_CASE(S8);