Changeset 797 for trunk/psLib/src/math/psStats.c

Timestamp:

May 27, 2004, 4:45:31 PM (22 years ago)

Author:

gusciora

Message:

More mods to get in sync with the latest SDRS.

File:

• trunk/psLib/src/math/psStats.c (modified) (23 diffs)

trunk/psLib/src/math/psStats.c

@@ -15 +15 @@
 
 #include "float.h"
-#define ROBUST_SIZE_THRESHOLD 10000   // Vectors that are large than this
+#define ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
 // will use robust statistical methods.
 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar
@@ -27 +27 @@
 #define false 0
 #define MYMAXFLOAT 1e99
+
+void p_psVectorRobustStats(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *stats);
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
@@ -40 +45 @@
     newStruct->sampleUQ = NAN;
     newStruct->sampleLQ = NAN;
-    newStruct->sampleLimit = NAN;
+    newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD;
     newStruct->robustMean = NAN;
     newStruct->robustMedian = NAN;
@@ -51 +56 @@
     newStruct->clippedMean = NAN;
     newStruct->clippedStdev = NAN;
-    newStruct->clipSigma = NAN;
-    newStruct->clipIter = -1;
+    newStruct->clipSigma = 3.0;
+    newStruct->clipIter = 3.0;
     newStruct->min = NAN;
     newStruct->max = NAN;
@@ -187 +192 @@
                    psVector *maskVector,
                    unsigned int maskVal,
-                   psStats *newStruct)
+                   psStats *stats)
 {
     int i = 0;
@@ -207 +212 @@
  ******************************************************************************
  *****************************************************************************/
+
+
+/******************************************************************************
+    ASSUMPTION: the mean is always calculated exactly.  Robust means are never
+    calculated in this routine.
+ *****************************************************************************/
 void p_psVectorSampleMean(const psVector *restrict myVector,
                           const psVector *restrict maskVector,
@@ -368 +379 @@
 
 /******************************************************************************
- This routine calculates the number of non-masked pixels in the vector. 
+    This routine calculates the number of non-masked pixels in the vector.
+    that fall within the min/max range, if given.
  *****************************************************************************/
 int p_psVectorNValues(const psVector *restrict myVector,
                       const psVector *restrict maskVector,
                       unsigned int maskVal,
-                      psStats *newStruct)
+                      psStats *stats)
 {
     int i = 0;
     int numData = 0;
-
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
-                numData++;
+    float rangeMin = 0.0;
+    float rangeMax = 0.0;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i]) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    numData++;
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    numData++;
+                }
             }
         }
     } else {
-        numData = myVector->n;
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    numData++;
+                }
+            }
+        } else {
+            numData = myVector->n;
+        }
     }
     return(numData);
@@ -395 +432 @@
                             const psVector *restrict maskVector,
                             unsigned int maskVal,
-                            psStats *newStruct)
+                            psStats *stats)
 {
     psVector *unsortedVector = NULL;
@@ -401 +438 @@
     int count = 0;
     int i = 0;
-    float median = 0.0;
     int nValues = 0;
-
-    nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct);
+    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_MEDIAN);
+        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
+
+        // Store the robust quartiles into the sample quartile members.
+        stats->sampleMedian = stats2->robustMedian;
+
+        // 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
+    // and are not maxed, IF the maskVector is not NULL>
+    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
+
+    // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues);
     unsortedVector->n = unsortedVector->nalloc;
@@ -410 +474 @@
     sortedVector->n = sortedVector->nalloc;
 
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
-                unsortedVector->vec.f[count++] = myVector->vec.f[i];
-
-            }
-        }
-        psSort(sortedVector, unsortedVector);
+    // Determine if we must only use data points within a min/max range.
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+
+        // Store all non-masked data points within the min/max range
+        // into the temporary vectors.
+        count = 0;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i]) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        }
     } else {
-        psSort(sortedVector, myVector);
-    }
-
+        // Store all non-masked data points into the temporary vectors.
+        count = 0;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->vec.f[i] = maskVector->vec.f[i];
+            }
+        }
+    }
+    // Sort the temporary vectors.
+    psSort(sortedVector, unsortedVector);
+
+    // Calculate the median exactly.
     if (0 == (nValues % 2)) {
-        median = 0.5 * (sortedVector->vec.f[(nValues/2)-1] +
-                        sortedVector->vec.f[nValues/2]);
+        stats->sampleMedian = 0.5 * (sortedVector->vec.f[(nValues/2)-1] +
+                                     sortedVector->vec.f[nValues/2]);
     } else {
-        median = sortedVector->vec.f[nValues/2];
-    }
-
+        stats->sampleMedian = sortedVector->vec.f[nValues/2];
+    }
+
+    // Free the temporary data structures.
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
-    newStruct->sampleMedian = median;
-}
-
+}
+
+/******************************************************************************
+    This routine smoothes the data in the input robustHistogram with a
+    Gaussian of width sigma.
+ *****************************************************************************/
 void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
                                   float sigma)
@@ -481 +580 @@
     p_psVectorSampleQuartiles()
  This procedure calculates the upper and/or lower quartiles of the
- data set.  It is assumed that the data set is small enough that we
- can sort all the data points and calculate the quartiles exactly.
+ data set.
+ 
  *****************************************************************************/
 void p_psVectorSampleQuartiles(const psVector *restrict myVector,
                                const psVector *restrict maskVector,
                                unsigned int maskVal,
-                               psStats *newStruct)
+                               psStats *stats)
 {
     psVector *unsortedVector = NULL;
@@ -495 +594 @@
     int i = 0;
     int nValues = 0;
-
-    // return is we have already calculated both quartile points.
-    if ((!isnan(newStruct->sampleLQ)) &&
-            (!isnan(newStruct->sampleUQ))) {
+    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;
     }
 
-    nValues = p_psVectorNValues(myVector, maskVector, maskVal, newStruct);
-
+    // Determine how many data points fit inside this min/max range
+    // and are not maxed, IF the maskVector is not NULL>
+    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
+
+    // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, nValues);
+    unsortedVector->n = unsortedVector->nalloc;
     sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, nValues);
-
-    count = 0;
-    if (maskVector != NULL) {
-        for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.ui8[i])) {
-                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
-            }
-        }
-        psSort(sortedVector, unsortedVector);
+    sortedVector->n = sortedVector->nalloc;
+
+    // Determine if we must only use data points within a min/max range.
+    if (stats->options & PS_STAT_USE_RANGE) {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        // Store all non-masked data points within the min/max range
+        // into the temporary vectors.
+        count = 0;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i]) &&
+                        (rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->vec.f[i]) &&
+                        (myVector->vec.f[i] <= rangeMax)) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        }
     } else {
-        psSort(sortedVector, myVector);
-    }
-
+        // Store all non-masked data points into the temporary vectors.
+        count = 0;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.ui8[i])) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                unsortedVector->vec.f[i] = maskVector->vec.f[i];
+            }
+        }
+    }
+
+    // Sort the temporary vectors.
+    psSort(sortedVector, unsortedVector);
+
+    // Calculate the quartile points exactly.
     ind = 3 * (nValues / 4);
-    newStruct->sampleUQ = sortedVector->vec.f[ind];
+    stats->sampleUQ = sortedVector->vec.f[ind];
     ind = (nValues / 4);
-    newStruct->sampleLQ = sortedVector->vec.f[ind];
-
+    stats->sampleLQ = sortedVector->vec.f[ind];
+
+    // Free the temporary data structures.
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
@@ -530 +683 @@
 
 /******************************************************************************
-    p_psVectorRobustStats(): this procedure calcualtes a variety of robust
+    p_psVectorRobustStats(): this procedure calculates a variety of robust
     stat measures:
- PS_STAT_ROBUST_MEAN
+        PS_STAT_ROBUST_MEAN
         PS_STAT_ROBUST_MEDIAN
         PS_STAT_ROBUST_MODE
         PS_STAT_ROBUST_STDEV
+        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
@@ -547 +701 @@
                            const psVector *restrict maskVector,
                            unsigned int maskVal,
-                           psStats *newStruct)
+                           psStats *stats)
 {
     psHistogram *robustHistogram = NULL;
@@ -565 +719 @@
     // is really required.
 
-    if (0.0 == newStruct->sampleUQ) {
+    if (0.0 == stats->sampleUQ) {
+
         // GUS: fix this
-        //        newStruct->options = newStruct->options | PS_STAT_SAMPLE_UQ;
+        //        stats->options = stats->options | PS_STAT_SAMPLE_UQ;
         p_psVectorSampleQuartiles(myVector,
                                   maskVector,
                                   maskVal,
-                                  newStruct);
-    }
-
-    if (isnan(newStruct->sampleLQ)) {
+                                  stats);
+    }
+
+    if (isnan(stats->sampleLQ)) {
         // GUS: fix this
-        //        newStruct->options = newStruct->options | PS_STAT_SAMPLE_LQ;
+        //        stats->options = stats->options | PS_STAT_SAMPLE_LQ;
         p_psVectorSampleQuartiles(myVector,
                                   maskVector,
                                   maskVal,
-                                  newStruct);
+                                  stats);
     }
 
     // Compute the initial bin size of the robust histogram.
-    sigmaE = (newStruct->sampleUQ - newStruct->sampleLQ) / 1.34f;
+    sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f;
     binSize = sigmaE / 10.0f;
 
     // Detemine minimum and maximum values in the data vector.
-    if (isnan(newStruct->min)) {
-        p_psVectorMin(myVector, maskVector, maskVal, newStruct);
-    }
-    if (isnan(newStruct->max)) {
-        p_psVectorMax(myVector, maskVector, maskVal, newStruct);
+    if (isnan(stats->min)) {
+        p_psVectorMin(myVector, maskVector, maskVal, stats);
+    }
+    if (isnan(stats->max)) {
+        p_psVectorMax(myVector, maskVector, maskVal, stats);
     }
 
     // Create the histogram structure.
-    robustHistogram = psHistogramAlloc(newStruct->min,
-                                       newStruct->max,
+    robustHistogram = psHistogramAlloc(stats->min,
+                                       stats->max,
                                        binSize);
     // Populate the histogram arrat.
@@ -610 +765 @@
     UQBinNum = -1;
     for (i=0;i<robustHistogram->nums->n;i++) {
-        if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleLQ) &&
-                (newStruct->sampleLQ <= robustHistogram->nums->vec.i32[i])) {
+        if ((robustHistogram->nums->vec.i32[i] <= stats->sampleLQ) &&
+                (stats->sampleLQ <= robustHistogram->nums->vec.i32[i])) {
             LQBinNum = i;
         }
 
-        if ((robustHistogram->nums->vec.i32[i] <= newStruct->sampleUQ) &&
-                (newStruct->sampleUQ <= robustHistogram->nums->vec.i32[i])) {
+        if ((robustHistogram->nums->vec.i32[i] <= stats->sampleUQ) &&
+                (stats->sampleUQ <= robustHistogram->nums->vec.i32[i])) {
             UQBinNum = i;
         }
@@ -631 +786 @@
     }
 
-    dL = (newStruct->robustUQ - newStruct->robustLQ) / 8.0;
+    dL = (stats->robustUQ - stats->robustLQ) / 8.0;
 
     // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
@@ -639 +794 @@
     // What is the mean_r?
 
-    if (newStruct->options & PS_STAT_ROBUST_MEAN) {
-        newStruct->robustMean = 0.0;
-    }
-
-    if  (newStruct->options & PS_STAT_ROBUST_MEDIAN) {
-        newStruct->robustMedian = 0.0;
-    }
-    if  (newStruct->options & PS_STAT_ROBUST_MODE) {
-        newStruct->robustMode = maxBinNum;
-    }
-    if  (newStruct->options & PS_STAT_ROBUST_STDEV) {
-        newStruct->robustStdev = 0.0;
-    }
-    if  (newStruct->options & PS_STAT_ROBUST_QUARTILE) {
-        newStruct->robustUQ = 0.0;
-        newStruct->robustLQ = 0.0;
+    if (stats->options & PS_STAT_ROBUST_MEAN) {
+        stats->robustMean = 0.0;
+    }
+
+    if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
+        stats->robustMedian = 0.0;
+    }
+    if  (stats->options & PS_STAT_ROBUST_MODE) {
+        stats->robustMode = maxBinNum;
+    }
+    if  (stats->options & PS_STAT_ROBUST_STDEV) {
+        stats->robustStdev = 0.0;
+    }
+    if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
+        stats->robustUQ = 0.0;
+        stats->robustLQ = 0.0;
     }
 }
@@ -661 +816 @@
     NOTE: This function assumes that p_psVectorMean() has already been called
     and the correct value is stored in stats->sampleMean.
+ 
+    NOTE: the mean is always calculated exactly.  Robust means are never
+    calculated in this routine.
  *****************************************************************************/
 void p_psVectorSampleStdev(const psVector *restrict myVector,
@@ -742 +900 @@
                             const psVector *restrict maskVector,
                             unsigned int maskVal,
-                            psStats *newStruct)
+                            psStats *stats)
 {
     int i = 0;
@@ -750 +908 @@
     psVector *tmpMask = NULL;
 
-    if (!((CLIPPED_NUM_ITER_LB <= newStruct->clipIter ) &&
-            (newStruct->clipIter <= CLIPPED_NUM_ITER_UB))) {
+    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
+            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
         psAbort(__func__, "Unallowed value for clipIter (%d).\n",
-                newStruct->clipIter);
-    }
-
-    if (!((CLIPPED_SIGMA_LB <= newStruct->clipSigma ) &&
-            (newStruct->clipSigma <= CLIPPED_SIGMA_UB))) {
+                stats->clipIter);
+    }
+
+    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
+            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
         psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
-                newStruct->clipSigma);
+                stats->clipSigma);
     }
 
@@ -770 +928 @@
 
     // 1. Compute the sample median.
-    p_psVectorSampleMedian(myVector, maskVector, maskVal, newStruct);
+    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
 
     // 2. Compute the sample standard deviation.
-    p_psVectorSampleStdev(myVector, maskVector, maskVal, newStruct);
+    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
 
     // 3. Use the sample median as the first estimator of the mean X.
-    clippedMean = newStruct->sampleMean;
+    clippedMean = stats->sampleMean;
 
     // 4. Use the sample stdev as the first estimator of the mean stdev.
-    clippedStdev = newStruct->sampleStdev;
+    clippedStdev = stats->sampleStdev;
 
     // 5. Repeat N times:
 
-    for (i=0;i<newStruct->clipIter;i++) {
+    for (i=0;i<stats->clipIter;i++) {
         for (j=0;j<myVector->n;j++) {
             // a) Exclude all values x_i for which |x_i - x| > K * stdev
             if ( fabs(myVector->vec.f[j] - clippedMean) >
-                    (newStruct->clipSigma * clippedStdev)) {
+                    (stats->clipSigma * clippedStdev)) {
                 tmpMask->vec.ui8[i] = 0xff;
             }
@@ -793 +951 @@
             // GUS: I should probably create a new struct here since the
             // following calls will overwrite any old values in sampleMean.
-            p_psVectorSampleMedian(myVector, tmpMask, maskVal, newStruct);
-            p_psVectorSampleStdev(myVector, tmpMask, maskVal, newStruct);
+            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
+            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
 
             // c) Use the new mean for x
-            clippedMean = newStruct->sampleMean;
+            clippedMean = stats->sampleMean;
 
             // d) Use the new stdev for stdev
-            clippedStdev = newStruct->sampleStdev;
+            clippedStdev = stats->sampleStdev;
         }
     }
 
     // 7. The last calcuated value of x is the cliped mean.
-    if (newStruct->options & PS_STAT_CLIPPED_MEAN) {
-        newStruct->clippedMean = clippedMean;
+    if (stats->options & PS_STAT_CLIPPED_MEAN) {
+        stats->clippedMean = clippedMean;
     }
 
     // 8. The last calcuated value of stdev is the cliped stdev.
-    if (newStruct->options & PS_STAT_CLIPPED_STDEV) {
-        newStruct->clippedStdev = clippedStdev;
+    if (stats->options & PS_STAT_CLIPPED_STDEV) {
+        stats->clippedStdev = clippedStdev;
     }