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

Timestamp:

Jun 6, 2004, 2:32:53 PM (22 years ago)

Author:

gusciora

Message:

Lots of bug fixes and such.

File:

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

trunk/psLib/src/math/psStats.c

@@ -85 +85 @@
     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 +98 @@
     }
     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 +109 @@
 }
 
-
+// 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 +118 @@
     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 +135 @@
 }
 
-void psHistogramFree(psHistogram *restrict myHist)
-{
+void psHistogramFree(psHistogram *myHist)
+{
+    psVectorFree(myHist->bounds);
+    psVectorFree(myHist->nums);
     psFree(myHist);
 }
@@ -150 +168 @@
 
     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 +225 @@
                                     binSize);
                     (out->nums->data.S32[binNum])++;
+
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
@@ -185 +240 @@
             }
         }
+
     }
     return(out);
@@ -471 +527 @@
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     unsortedVector->n = unsortedVector->nalloc;
+
     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     sortedVector->n = sortedVector->nalloc;
@@ -494 +551 @@
                 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 +561 @@
             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 +653 @@
     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 +699 @@
                         (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 +706 @@
                 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 +716 @@
             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 +738 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
+    // GUS: This is the
 }
 
@@ -691 +751 @@
         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 +779 @@
     // 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 +788 @@
     }
 
-    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 +803 @@
     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 +863 @@
         stats->robustLQ = 0.0;
     }
+    stats->robustNfit = 0.0;
+    stats->robustN50 = 0.0;
+    psHistogramFree(robustHistogram);
 }
 
@@ -836 +891 @@
 
     if (0 != isnan(stats->sampleMean)) {
-        psAbort(__func__, "stats->sampleMean == NAN");
+        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     }
     mean = stats->sampleMean;
@@ -920 +975 @@
     }
 
-    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 +1002 @@
 
     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 +1056 @@
     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 +1082 @@
 
     // ************************************************************************
+    // 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);
     }