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

Timestamp:

Aug 4, 2004, 1:37:39 PM (22 years ago)

Author:

desonia

Message:

found the server astyle upgrade was faulty, so the format was reset.

File:

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

trunk/psLib/src/math/psStats.c

@@ -36 +36 @@
                             unsigned int maskVal,
                             psStats *stats );
-                            
+
 /** Preprocessor macro to generate error on an incorrect type */
 #define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \
 if (NAME->type.type != TYPE) { \
-        psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \
-    }
-    
+    psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \
+}
+
 /** Preprocessor macro to generate error on a NULL vector */
 #define PS_CHECK_NULL_VECTOR(NAME) \
 if (NAME == NULL || NAME->data.V == NULL) { \
-        psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
-    }
-    
+    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
+}
+
 /** Preprocessor macro to generate error for zero length vector */
 #define PS_CHECK_EMPTY_VECTOR(NAME) \
 if (NAME->n < 1) { \
-        psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
-    }
-    
+    psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+}
+
 /** Preprocessor macro to generate error on differing size vectors */
 #define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \
 if (VEC1->n != VEC2->n) { \
-        psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
-    }
-    
+    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
+}
+
 #define PS_PRINT_VECTOR(NAME) \
 for (int my_i=0;my_i<NAME->n;my_i++) { \
-        printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \
-    } \
+    printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \
+} \
 printf("\n"); \
 
@@ -92 +92 @@
     switch ( stats->options &
              ~ ( PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE ) ) {
-            case PS_STAT_SAMPLE_MEAN:
-            *value = stats->sampleMean;
-            return true;
-            
-            case PS_STAT_SAMPLE_MEDIAN:
-            *value = stats->sampleMedian;
-            return true;
-            
-            case PS_STAT_SAMPLE_STDEV:
-            *value = stats->sampleStdev;
-            return true;
-            
-            case PS_STAT_ROBUST_MEAN:
-            *value = stats->robustMean;
-            return true;
-            
-            case PS_STAT_ROBUST_MEDIAN:
-            *value = stats->robustMedian;
-            return true;
-            
-            case PS_STAT_ROBUST_MODE:
-            *value = stats->robustMode;
-            return true;
-            
-            case PS_STAT_ROBUST_STDEV:
-            *value = stats->robustStdev;
-            return true;
-            
-            case PS_STAT_CLIPPED_MEAN:
-            *value = stats->clippedMean;
-            return true;
-            
-            case PS_STAT_CLIPPED_STDEV:
-            *value = stats->clippedStdev;
-            return true;
-            
-            case PS_STAT_MAX:
-            *value = stats->max;
-            return true;
-            
-            case PS_STAT_MIN:
-            *value = stats->min;
-            return true;
-            
-            default:
-            return false;
-        }
+    case PS_STAT_SAMPLE_MEAN:
+        *value = stats->sampleMean;
+        return true;
+
+    case PS_STAT_SAMPLE_MEDIAN:
+        *value = stats->sampleMedian;
+        return true;
+
+    case PS_STAT_SAMPLE_STDEV:
+        *value = stats->sampleStdev;
+        return true;
+
+    case PS_STAT_ROBUST_MEAN:
+        *value = stats->robustMean;
+        return true;
+
+    case PS_STAT_ROBUST_MEDIAN:
+        *value = stats->robustMedian;
+        return true;
+
+    case PS_STAT_ROBUST_MODE:
+        *value = stats->robustMode;
+        return true;
+
+    case PS_STAT_ROBUST_STDEV:
+        *value = stats->robustStdev;
+        return true;
+
+    case PS_STAT_CLIPPED_MEAN:
+        *value = stats->clippedMean;
+        return true;
+
+    case PS_STAT_CLIPPED_STDEV:
+        *value = stats->clippedStdev;
+        return true;
+
+    case PS_STAT_MAX:
+        *value = stats->max;
+        return true;
+
+    case PS_STAT_MIN:
+        *value = stats->min;
+        return true;
+
+    default:
+        return false;
+    }
 }
 
@@ -153 +153 @@
 {
     int i = 0;                                  // Loop index variable.
-    
+
     for ( i = 0;i < myVector->n;i++ ) {
-            if ( maskVector != NULL )
-                printf( "Element %d is %f (mask is %d)\n", i, myVector->data.F32[ i ], maskVector->data.U8[ i ] );
-            else
-                printf( "Element %d is %f\n", i, myVector->data.F32[ i ] );
-        }
+        if ( maskVector != NULL )
+            printf( "Element %d is %f (mask is %d)\n", i, myVector->data.F32[ i ], maskVector->data.U8[ i ] );
+        else
+            printf( "Element %d is %f\n", i, myVector->data.F32[ i ] );
+    }
 }
 
@@ -204 +204 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
     // loop.
     if ( stats->options & PS_STAT_USE_RANGE ) {
-            rangeMin = stats->min;
-            rangeMax = stats->max;
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            // Check if the data is with the specified range
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    mean += myVector->data.F32[ i ];
-                                    count++;
-                                }
-                        }
-                    mean /= ( float ) count;
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    mean += myVector->data.F32[ i ];
-                                    count++;
-                                }
-                        }
-                    mean /= ( float ) count;
-                }
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                // Check if the data is with the specified range
+                if ( !( maskVal & maskVector->data.U8[ i ] ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    mean += myVector->data.F32[ i ];
+                    count++;
+                }
+            }
+            mean /= ( float ) count;
         } else {
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
-                                    mean += myVector->data.F32[ i ];
-                                    count++;
-                                }
-                        }
-                    mean /= ( float ) count;
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            mean += myVector->data.F32[ i ];
-                        }
-                    mean /= ( float ) myVector->n;
-                }
-        }
-        
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    mean += myVector->data.F32[ i ];
+                    count++;
+                }
+            }
+            mean /= ( float ) count;
+        }
+    } else {
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    mean += myVector->data.F32[ i ];
+                    count++;
+                }
+            }
+            mean /= ( float ) count;
+        } else {
+            for ( i = 0;i < myVector->n;i++ ) {
+                mean += myVector->data.F32[ i ];
+            }
+            mean /= ( float ) myVector->n;
+        }
+    }
+
     stats->sampleMean = mean;
 }
@@ -271 +271 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
     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->data.U8[ i ] ) ) {
-                                    if ( ( myVector->data.F32[ i ] > max ) &&
-                                            ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                            ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                            max = myVector->data.F32[ i ];
-                                        }
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( myVector->data.F32[ i ] > max ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    max = myVector->data.F32[ i ];
-                                }
-                        }
-                }
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    if ( ( myVector->data.F32[ i ] > max ) &&
+                            ( rangeMin <= myVector->data.F32[ i ] ) &&
+                            ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                        max = myVector->data.F32[ i ];
+                    }
+                }
+            }
         } else {
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
-                                    if ( myVector->data.F32[ i ] > max ) {
-                                            max = myVector->data.F32[ i ];
-                                        }
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( myVector->data.F32[ i ] > max ) {
-                                    max = myVector->data.F32[ i ];
-                                }
-                        }
-                }
-        }
-        
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( myVector->data.F32[ i ] > max ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    max = myVector->data.F32[ i ];
+                }
+            }
+        }
+    } else {
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    if ( myVector->data.F32[ i ] > max ) {
+                        max = myVector->data.F32[ i ];
+                    }
+                }
+            }
+        } else {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( myVector->data.F32[ i ] > max ) {
+                    max = myVector->data.F32[ i ];
+                }
+            }
+        }
+    }
+
     stats->max = max;
 }
@@ -336 +336 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     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->data.U8[ i ] ) ) {
-                                    if ( ( myVector->data.F32[ i ] < min ) &&
-                                            ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                            ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                            min = myVector->data.F32[ i ];
-                                        }
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( myVector->data.F32[ i ] < min ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    min = myVector->data.F32[ i ];
-                                }
-                        }
-                }
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    if ( ( myVector->data.F32[ i ] < min ) &&
+                            ( rangeMin <= myVector->data.F32[ i ] ) &&
+                            ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                        min = myVector->data.F32[ i ];
+                    }
+                }
+            }
         } else {
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
-                                    if ( myVector->data.F32[ i ] < min ) {
-                                            min = myVector->data.F32[ i ];
-                                        }
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( myVector->data.F32[ i ] < min ) {
-                                    min = myVector->data.F32[ i ];
-                                }
-                        }
-                }
-        }
-        
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( myVector->data.F32[ i ] < min ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    min = myVector->data.F32[ i ];
+                }
+            }
+        }
+    } else {
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    if ( myVector->data.F32[ i ] < min ) {
+                        min = myVector->data.F32[ i ];
+                    }
+                }
+            }
+        } else {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( myVector->data.F32[ i ] < min ) {
+                    min = myVector->data.F32[ i ];
+                }
+            }
+        }
+    }
+
     stats->min = min;
 }
@@ -401 +401 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     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->data.U8[ i ] ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    numData++;
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    numData++;
-                                }
-                        }
-                }
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    numData++;
+                }
+            }
         } else {
-            rangeMin = stats->min;
-            rangeMax = stats->max;
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
-                                    numData++;
-                                }
-                        }
-                } else {
-                    numData = myVector->n;
-                }
-        }
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    numData++;
+                }
+            }
+        }
+    } else {
+        rangeMin = stats->min;
+        rangeMax = stats->max;
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    numData++;
+                }
+            }
+        } else {
+            numData = myVector->n;
+        }
+    }
     return ( numData );
 }
@@ -462 +462 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
-    
+
+
     // Determine if the number of data points exceed a threshold which will
     // cause to generate robust stats, as opposed to exact stats.
@@ -488 +488 @@
         }
     */
-    
+
     // Determine how many data points fit inside this min/max range
     // and are not masked, IF the maskVector is not NULL>
     nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats );
-    
+
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     unsortedVector->n = unsortedVector->nalloc;
-    
+
     sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     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->data.U8[ i ] ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    unsortedVector->data.F32[ count++ ] = maskVector->data.F32[ i ];
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
-                                }
-                        }
-                }
+        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->data.U8[ i ] ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    unsortedVector->data.F32[ count++ ] = maskVector->data.F32[ i ];
+                }
+            }
         } else {
-            // 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->data.U8[ 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 ];
-                        }
-                }
-        }
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
+                }
+            }
+        }
+    } else {
+        // 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->data.U8[ 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 ];
+            }
+        }
+    }
     // Sort the temporary vectors.
     psVectorSort( sortedVector, unsortedVector );
-    
+
     // Calculate the median exactly.
     // XXX: Is this the correct action?
     if ( 0 == ( nValues % 2 ) ) {
-            stats->sampleMedian = 0.5 * ( sortedVector->data.F32[ ( nValues / 2 ) - 1 ] +
-                                          sortedVector->data.F32[ nValues / 2 ] );
-        } else {
-            stats->sampleMedian = sortedVector->data.F32[ nValues / 2 ];
-        }
-        
+        stats->sampleMedian = 0.5 * ( sortedVector->data.F32[ ( nValues / 2 ) - 1 ] +
+                                      sortedVector->data.F32[ nValues / 2 ] );
+    } else {
+        stats->sampleMedian = sortedVector->data.F32[ nValues / 2 ];
+    }
+
     // Free the temporary data structures.
     psFree( unsortedVector );
@@ -577 +577 @@
     float gaussianCoefs[ 1 + ( 2 * GAUSS_WIDTH ) ]; // The Gaussian Coefficients
     psVector *smooth = psVectorAlloc( robustHistogram->nums->n, PS_TYPE_F32 );
-    
+
     for ( i = 0;i < ( 1 + ( 2 * GAUSS_WIDTH ) );i++ ) {
-            if ( fabs( sigma ) >= FLT_EPSILON ) {
-                    // If sigma does not equal zero, then we use Gaussian smoothing.
-                    #ifdef  DARWIN
-                    denom = ( float ) sqrt( 2.0 * M_PI * sigma * sigma );
-                    #else
-                    
-                    denom = sqrtf( 2.0 * M_PI * sigma * sigma );
-                    #endif
-                    
-                    expo = - ( float ) ( ( i - GAUSS_WIDTH ) * ( i - GAUSS_WIDTH ) );
-                    expo /= ( 2.0 * sigma * sigma );
-                    gaussianCoefs[ i ] = exp( expo / denom );
-                    
-                    // NOTE: Gaussian smoothing just isn't working with low sigma
-                    // values.  The problem is that the Gaussian coefficients are
-                    // all zero, except for the middle coefficient, which is exactly
-                    // one.  Therefore, I'm using boxcar smoothing.
-                    gaussianCoefs[ i ] = 1.0 / ( 1.0 + ( 2.0 * ( float ) GAUSS_WIDTH ) );
-                    //            printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-                } else {
-                    /* If sigma equals zero (all pixels have the same value)
-                     * the above code will divide by zero.  Therefore, we don't need
-                     * to smooth the data.
-                     */
-                    for ( i = 0;i < robustHistogram->nums->n;i++ ) {
-                            smooth->data.F32[ i ] = ( float ) robustHistogram->nums->data.S32[ i ];
-                        }
-                    return ( smooth );
-                }
-        }
-        
+        if ( fabs( sigma ) >= FLT_EPSILON ) {
+            // If sigma does not equal zero, then we use Gaussian smoothing.
+            #ifdef  DARWIN
+            denom = ( float ) sqrt( 2.0 * M_PI * sigma * sigma );
+            #else
+
+            denom = sqrtf( 2.0 * M_PI * sigma * sigma );
+            #endif
+
+            expo = - ( float ) ( ( i - GAUSS_WIDTH ) * ( i - GAUSS_WIDTH ) );
+            expo /= ( 2.0 * sigma * sigma );
+            gaussianCoefs[ i ] = exp( expo / denom );
+
+            // NOTE: Gaussian smoothing just isn't working with low sigma
+            // values.  The problem is that the Gaussian coefficients are
+            // all zero, except for the middle coefficient, which is exactly
+            // one.  Therefore, I'm using boxcar smoothing.
+            gaussianCoefs[ i ] = 1.0 / ( 1.0 + ( 2.0 * ( float ) GAUSS_WIDTH ) );
+            //            printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
+        } else {
+            /* If sigma equals zero (all pixels have the same value)
+             * the above code will divide by zero.  Therefore, we don't need
+             * to smooth the data.
+             */
+            for ( i = 0;i < robustHistogram->nums->n;i++ ) {
+                smooth->data.F32[ i ] = ( float ) robustHistogram->nums->data.S32[ i ];
+            }
+            return ( smooth );
+        }
+    }
+
     // Perform the actual smoothing.
     for ( i = 0;i < robustHistogram->nums->n;i++ ) {
-            smooth->data.F32[ i ] = 0.0;
-            for ( j = -GAUSS_WIDTH;j <= + GAUSS_WIDTH;j++ ) {
-                    if ( ( ( j + i ) >= 0 ) && ( ( j + i ) < smooth->n ) ) {
-                            smooth->data.F32[ i ] += ( gaussianCoefs[ j + GAUSS_WIDTH ] *
-                                                       ( float ) robustHistogram->nums->data.S32[ j + i ] );
-                        }
-                }
-        }
+        smooth->data.F32[ i ] = 0.0;
+        for ( j = -GAUSS_WIDTH;j <= + GAUSS_WIDTH;j++ ) {
+            if ( ( ( j + i ) >= 0 ) && ( ( j + i ) < smooth->n ) ) {
+                smooth->data.F32[ i ] += ( gaussianCoefs[ j + GAUSS_WIDTH ] *
+                                           ( float ) robustHistogram->nums->data.S32[ j + i ] );
+            }
+        }
+    }
     return ( smooth );
 }
@@ -646 +646 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     // 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( nValues, PS_TYPE_F32 );
@@ -656 +656 @@
     sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     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->data.U8[ i ] ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
-                                }
-                        }
-                }
+        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->data.U8[ i ] ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
+                }
+            }
         } else {
-            // 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->data.U8[ 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 ];
-                        }
-                }
-        }
-        
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
+                }
+            }
+        }
+    } else {
+        // 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->data.U8[ 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 ];
+            }
+        }
+    }
+
     // Sort the temporary vectors.
     psVectorSort( sortedVector, unsortedVector );
-    
+
     // Calculate the quartile points exactly.
     stats->sampleUQ = sortedVector->data.F32[ 3 * ( nValues / 4 ) ];
     stats->sampleLQ = sortedVector->data.F32[ nValues / 4 ];
-    
+
     // Free the temporary data structures.
     psFree( unsortedVector );
@@ -735 +735 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    
+
     // This procedure requires the mean.  If it has not been already
     // calculated, then call p_psVectorSampleMean()
     if ( 0 != isnan( stats->sampleMean ) ) {
-            p_psVectorSampleMean( myVector, maskVector, maskVal, stats );
-        }
+        p_psVectorSampleMean( myVector, maskVector, maskVal, stats );
+    }
     mean = stats->sampleMean;
-    
+
     if ( stats->options & PS_STAT_USE_RANGE ) {
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) &&
-                                    ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    diff = myVector->data.F32[ i ] - mean;
-                                    sumSquares += ( diff * diff );
-                                    sumDiffs += diff;
-                                    countInt++;
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
-                                    ( myVector->data.F32[ i ] <= rangeMax ) ) {
-                                    diff = myVector->data.F32[ i ] - mean;
-                                    sumSquares += ( diff * diff );
-                                    sumDiffs += diff;
-                                    countInt++;
-                                }
-                        }
-                    countInt = myVector->n;
-                }
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) &&
+                        ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    diff = myVector->data.F32[ i ] - mean;
+                    sumSquares += ( diff * diff );
+                    sumDiffs += diff;
+                    countInt++;
+                }
+            }
         } else {
-            if ( maskVector != NULL ) {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
-                                    diff = myVector->data.F32[ i ] - mean;
-                                    sumSquares += ( diff * diff );
-                                    sumDiffs += diff;
-                                    countInt++;
-                                }
-                        }
-                } else {
-                    for ( i = 0;i < myVector->n;i++ ) {
-                            diff = myVector->data.F32[ i ] - mean;
-                            sumSquares += ( diff * diff );
-                            sumDiffs += diff;
-                            countInt++;
-                        }
-                    countInt = myVector->n;
-                }
-        }
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
+                        ( myVector->data.F32[ i ] <= rangeMax ) ) {
+                    diff = myVector->data.F32[ i ] - mean;
+                    sumSquares += ( diff * diff );
+                    sumDiffs += diff;
+                    countInt++;
+                }
+            }
+            countInt = myVector->n;
+        }
+    } else {
+        if ( maskVector != NULL ) {
+            for ( i = 0;i < myVector->n;i++ ) {
+                if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
+                    diff = myVector->data.F32[ i ] - mean;
+                    sumSquares += ( diff * diff );
+                    sumDiffs += diff;
+                    countInt++;
+                }
+            }
+        } else {
+            for ( i = 0;i < myVector->n;i++ ) {
+                diff = myVector->data.F32[ i ] - mean;
+                sumSquares += ( diff * diff );
+                sumDiffs += diff;
+                countInt++;
+            }
+            countInt = myVector->n;
+        }
+    }
     countFloat = ( float ) countInt;
-    
+
     #ifdef DARWIN
-    
+
     stats->sampleStdev = ( float ) sqrt( ( sumSquares - ( sumDiffs *
                                            sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     #else
-    
+
     stats->sampleStdev = sqrtf( ( sumSquares - ( sumDiffs *
                                   sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
@@ -824 +824 @@
     float oldStanStdev = 0.0;                   // Temporary variable
     psVector *tmpMask = NULL;                   // Temporary vector
-    
+
     // Endure that stats->clipIter is within the proper range.
     if ( !( ( CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
             ( stats->clipIter <= CLIPPED_NUM_ITER_UB ) ) ) {
-            psAbort( __func__, "Unallowed value for clipIter (%d).\n",
-                     stats->clipIter );
-        }
-        
+        psAbort( __func__, "Unallowed value for clipIter (%d).\n",
+                 stats->clipIter );
+    }
+
     // Endure that stats->clipSigma is within the proper range.
     if ( !( ( CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
             ( stats->clipSigma <= CLIPPED_SIGMA_UB ) ) ) {
-            psAbort( __func__, "Unallowed value for clipSigma (%f).\n",
-                     stats->clipSigma );
-        }
-        
+        psAbort( __func__, "Unallowed value for clipSigma (%f).\n",
+                 stats->clipSigma );
+    }
+
     // We allocate a temporary mask vector since during the iterative
     // steps that follow, we will be masking off additional data points.
@@ -844 +844 @@
     tmpMask = psVectorAlloc( myVector->n, PS_TYPE_U8 );
     tmpMask->n = myVector->n;
-    
+
     // If we were called with a mask vector, then initialize the temporary
     // mask vector with those values.
     if ( maskVector != NULL ) {
-            for ( i = 0;i < tmpMask->n;i++ ) {
-                    tmpMask->data.U8[ i ] = maskVector->data.U8[ i ];
-                }
-        }
-        
+        for ( i = 0;i < tmpMask->n;i++ ) {
+            tmpMask->data.U8[ i ] = maskVector->data.U8[ i ];
+        }
+    }
+
     // 1. Compute the sample median.
     // NOTE: This seems odd.  Verify with IfA that we want to calculate the
     // median here, not the mean.
     p_psVectorSampleMedian( myVector, maskVector, maskVal, stats );
-    
+
     // 2. Compute the sample standard deviation.
     p_psVectorSampleStdev( myVector, maskVector, maskVal, stats );
-    
+
     // 3. Use the sample median as the first estimator of the mean X.
     clippedMean = stats->sampleMean;
-    
+
     // 4. Use the sample stdev as the first estimator of the mean stdev.
     clippedStdev = stats->sampleStdev;
-    
+
     // Must save the old sampleMean and sampleStdev since the following code
     // block overwrites them.
     oldStanMean = stats->sampleMean;
     oldStanStdev = stats->sampleStdev;
-    
+
     // 5. Repeat N times:
     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->data.F32[ j ] - clippedMean ) >
-                            ( stats->clipSigma * clippedStdev ) ) {
-                            tmpMask->data.U8[ i ] = 0xff;
-                        }
-                    // b) compute new mean and stdev
-                    p_psVectorSampleMedian( myVector, tmpMask, maskVal, stats );
-                    p_psVectorSampleStdev( myVector, tmpMask, maskVal, stats );
-                    
-                    // c) Use the new mean for x
-                    clippedMean = stats->sampleMean;
-                    
-                    // d) Use the new stdev for stdev
-                    clippedStdev = stats->sampleStdev;
-                }
-                
-        }
+        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 ) >
+                    ( stats->clipSigma * clippedStdev ) ) {
+                tmpMask->data.U8[ i ] = 0xff;
+            }
+            // b) compute new mean and stdev
+            p_psVectorSampleMedian( myVector, tmpMask, maskVal, stats );
+            p_psVectorSampleStdev( myVector, tmpMask, maskVal, stats );
+
+            // c) Use the new mean for x
+            clippedMean = stats->sampleMean;
+
+            // d) Use the new stdev for stdev
+            clippedStdev = stats->sampleStdev;
+        }
+
+    }
     stats->sampleMean = oldStanMean;
     stats->sampleStdev = oldStanStdev;
-    
+
     // 7. The last calcuated value of x is the cliped mean.
     if ( stats->options & PS_STAT_CLIPPED_MEAN ) {
-            stats->clippedMean = clippedMean;
-        }
-        
+        stats->clippedMean = clippedMean;
+    }
+
     // 8. The last calcuated value of stdev is the cliped stdev.
     if ( stats->options & PS_STAT_CLIPPED_STDEV ) {
-            stats->clippedStdev = clippedStdev;
-        }
-        
+        stats->clippedStdev = clippedStdev;
+    }
+
     psFree( tmpMask );
 }
@@ -918 +918 @@
     float range = 0.0;
     int i = 0;
-    
+
     for ( i = 0;i < myData->n;i++ ) {
-            if ( myData->data.F32[ i ] < min ) {
-                    min = myData->data.F32[ i ];
-                }
-            if ( myData->data.F32[ i ] > max ) {
-                    max = myData->data.F32[ i ];
-                }
-        }
-        
+        if ( myData->data.F32[ i ] < min ) {
+            min = myData->data.F32[ i ];
+        }
+        if ( myData->data.F32[ i ] > max ) {
+            max = myData->data.F32[ i ];
+        }
+    }
+
     range = max - min;
     for ( i = 0;i < myData->n;i++ ) {
-            myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range;
-        }
+        myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range;
+    }
 }
 
@@ -948 +948 @@
     float tmp = exp( -( ( x - mean ) * ( x - mean ) ) / ( 2.0 * stdev * stdev ) );
     tmp /= ( ( float ) sqrt( 2.0 * M_PI * ( stdev * stdev ) ) );
-    
+
     //    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     return ( tmp );
@@ -965 +965 @@
     float stdev = myParams->data.F32[ 1 ];
     float tmp = 0.0;
-    
+
     if ( whichParam == 0 ) {
-            // Return the derivative w.r.t. the mean.
-            tmp = ( x - mean ) * p_psGaussian( myData, myParams );
-            tmp /= ( stdev * stdev );
-        } else if ( whichParam == 1 ) {
+        // Return the derivative w.r.t. the mean.
+        tmp = ( x - mean ) * p_psGaussian( myData, myParams );
+        tmp /= ( stdev * stdev );
+    } else
+        if ( whichParam == 1 ) {
             // Return the derivative w.r.t. the stdev.
             tmp = ( x - mean ) * ( x - mean ) * p_psGaussian( myData, myParams );
@@ -976 +977 @@
         }
     printf( "p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp );
-    
+
     return ( tmp );
 }
@@ -994 +995 @@
     float C = myParams->data.F32[ 2 ];
     float tmp = 0.0;
-    
+
     tmp = ( A * x * x ) + ( B * x ) + C;
     return ( tmp );
@@ -1009 +1010 @@
     float x = myCoords->data.F32[ 0 ];
     float tmp = 0.0;
-    
+
     if ( whichParamDeriv == 0 ) {
-            tmp = x * x;
-        } else if ( whichParamDeriv == 1 ) {
+        tmp = x * x;
+    } else
+        if ( whichParamDeriv == 1 ) {
             tmp = x;
-        } else if ( whichParamDeriv == 2 ) {
-            tmp = 1.0;
-        }
-        
+        } else
+            if ( whichParamDeriv == 2 ) {
+                tmp = 1.0;
+            }
+
     return ( tmp );
 }
@@ -1039 +1042 @@
     float oldMidpoint = 1.0;
     float f = 0.0;
-    
+
     //  printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
-    
+
     while ( numIterations < MAX_ITERATIONS ) {
-            midpoint = ( rangeHigh + rangeLow ) / 2.0;
-            if ( fabs( midpoint - oldMidpoint ) <= FLT_EPSILON ) {
-                    return ( midpoint );
-                }
-            oldMidpoint = midpoint;
-            
-            f = psPolynomial1DEval( midpoint, myPoly );
-            //        printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
-            if ( fabs( f - getThisValue ) <= FLT_EPSILON ) {
-                    return ( midpoint );
-                }
-                
-            if ( f > getThisValue ) {
-                    rangeHigh = midpoint;
-                } else {
-                    rangeLow = midpoint;
-                }
-            numIterations++;
-        }
+        midpoint = ( rangeHigh + rangeLow ) / 2.0;
+        if ( fabs( midpoint - oldMidpoint ) <= FLT_EPSILON ) {
+            return ( midpoint );
+        }
+        oldMidpoint = midpoint;
+
+        f = psPolynomial1DEval( midpoint, myPoly );
+        //        printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
+        if ( fabs( f - getThisValue ) <= FLT_EPSILON ) {
+            return ( midpoint );
+        }
+
+        if ( f > getThisValue ) {
+            rangeHigh = midpoint;
+        } else {
+            rangeLow = midpoint;
+        }
+        numIterations++;
+    }
     return ( midpoint );
 }
@@ -1082 +1085 @@
     psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 );
     psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 );
-    
+
     if ( ( binNum > 0 ) &&
             ( binNum < ( histogram->nums->n + 1 ) ) ) {
-            x->data.F64[ 0 ] = ( double ) 0.5 *
-                               ( histogram->bounds->data.F32[ binNum - 1 ] +
-                                 histogram->bounds->data.F32[ binNum ] );
-            x->data.F64[ 1 ] = ( double ) 0.5 *
-                               ( histogram->bounds->data.F32[ binNum ] +
-                                 histogram->bounds->data.F32[ binNum + 1 ] );
-            x->data.F64[ 2 ] = ( double ) 0.5 *
-                               ( histogram->bounds->data.F32[ binNum + 1 ] +
-                                 histogram->bounds->data.F32[ binNum + 2 ] );
-                                 
-            y->data.F64[ 0 ] = cumulativeSums->data.F32[ binNum - 1 ];
-            y->data.F64[ 1 ] = cumulativeSums->data.F32[ binNum ];
-            y->data.F64[ 2 ] = cumulativeSums->data.F32[ binNum + 1 ];
-            
-            if ( !( ( y->data.F64[ 0 ] <= fitFloat ) &&
-                    ( fitFloat <= y->data.F64[ 2 ] ) ) ) {
-                    psAbort( __func__, "p_psVectorRobustStats(0): midpoint not within y-range\n" );
-                }
-                
-            yErr->data.F64[ 0 ] = 1.0;
-            yErr->data.F64[ 1 ] = 1.0;
-            yErr->data.F64[ 2 ] = 1.0;
-            
-            myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
-            return ( p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], x->data.F64[ 2 ],
-                                       fitFloat ) );
-        } else {
-            return ( 0.5 * ( histogram->bounds->data.F32[ binNum + 1 ] +
-                             histogram->bounds->data.F32[ binNum ] ) );
-        }
-        
+        x->data.F64[ 0 ] = ( double ) 0.5 *
+                           ( histogram->bounds->data.F32[ binNum - 1 ] +
+                             histogram->bounds->data.F32[ binNum ] );
+        x->data.F64[ 1 ] = ( double ) 0.5 *
+                           ( histogram->bounds->data.F32[ binNum ] +
+                             histogram->bounds->data.F32[ binNum + 1 ] );
+        x->data.F64[ 2 ] = ( double ) 0.5 *
+                           ( histogram->bounds->data.F32[ binNum + 1 ] +
+                             histogram->bounds->data.F32[ binNum + 2 ] );
+
+        y->data.F64[ 0 ] = cumulativeSums->data.F32[ binNum - 1 ];
+        y->data.F64[ 1 ] = cumulativeSums->data.F32[ binNum ];
+        y->data.F64[ 2 ] = cumulativeSums->data.F32[ binNum + 1 ];
+
+        if ( !( ( y->data.F64[ 0 ] <= fitFloat ) &&
+                ( fitFloat <= y->data.F64[ 2 ] ) ) ) {
+            psAbort( __func__, "p_psVectorRobustStats(0): midpoint not within y-range\n" );
+        }
+
+        yErr->data.F64[ 0 ] = 1.0;
+        yErr->data.F64[ 1 ] = 1.0;
+        yErr->data.F64[ 2 ] = 1.0;
+
+        myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
+        return ( p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], x->data.F64[ 2 ],
+                                   fitFloat ) );
+    } else {
+        return ( 0.5 * ( histogram->bounds->data.F32[ binNum + 1 ] +
+                         histogram->bounds->data.F32[ binNum ] ) );
+    }
+
     psFree( x );
     psFree( y );
@@ -1187 +1190 @@
     float sumNfit = 0.0;
     float cumulativeMedian = 0.0;
-    
+
     // Compute the initial bin size of the robust histogram.  This is done
     // by computing the clipped standard deviation of the vector, and dividing
@@ -1193 +1196 @@
     p_psVectorClippedStats( myVector, maskVector, maskVal, tmpStats );
     binSize = tmpStats->clippedStdev / 10.0f;
-    
+
     // If stats->clippedStdev == 0.0, then all data elements have the same
     // value.  Therefore, we can set the appropiate results and return.
     if ( fabs( binSize ) <= FLT_EPSILON ) {
-            if ( stats->options & PS_STAT_ROBUST_MEAN ) {
-                    stats->robustMean = stats->clippedMean;
-                }
-            if ( stats->options & PS_STAT_ROBUST_MEDIAN ) {
-                    stats->robustMedian = stats->clippedMean;
-                }
-            if ( stats->options & PS_STAT_ROBUST_MODE ) {
-                    stats->robustMode = stats->clippedMean;
-                }
-            if ( stats->options & PS_STAT_ROBUST_STDEV ) {
-                    stats->robustStdev = 0.0;
-                }
-            if ( stats->options & PS_STAT_ROBUST_QUARTILE ) {
-                    stats->robustUQ = stats->clippedMean;
-                    stats->robustLQ = stats->clippedMean;
-                }
-            // XXX: Set these to the number of unmasked data points?
-            stats->robustNfit = 0.0;
-            stats->robustN50 = 0.0;
-            psFree( tmpStats );
-            return ;
-        }
-        
+        if ( stats->options & PS_STAT_ROBUST_MEAN ) {
+            stats->robustMean = stats->clippedMean;
+        }
+        if ( stats->options & PS_STAT_ROBUST_MEDIAN ) {
+            stats->robustMedian = stats->clippedMean;
+        }
+        if ( stats->options & PS_STAT_ROBUST_MODE ) {
+            stats->robustMode = stats->clippedMean;
+        }
+        if ( stats->options & PS_STAT_ROBUST_STDEV ) {
+            stats->robustStdev = 0.0;
+        }
+        if ( stats->options & PS_STAT_ROBUST_QUARTILE ) {
+            stats->robustUQ = stats->clippedMean;
+            stats->robustLQ = stats->clippedMean;
+        }
+        // XXX: Set these to the number of unmasked data points?
+        stats->robustNfit = 0.0;
+        stats->robustN50 = 0.0;
+        psFree( tmpStats );
+        return ;
+    }
+
     // Determine minimum and maximum values in the data vector.
     if ( isnan( stats->min ) ) {
-            p_psVectorMin( myVector, maskVector, maskVal, stats );
-        }
+        p_psVectorMin( myVector, maskVector, maskVal, stats );
+    }
     if ( isnan( stats->max ) ) {
-            p_psVectorMax( myVector, maskVector, maskVal, stats );
-        }
-        
+        p_psVectorMax( myVector, maskVector, maskVal, stats );
+    }
+
     // Create the histogram structure.  NOTE: we can not specify the bin size
     // precisely since the argument to psHistogramAlloc() is the number of
@@ -1236 +1239 @@
                                         stats->max,
                                         numBins );
-                                        
+
     // Populate the histogram array.
     psVectorHistogram( robustHistogram, myVector, maskVector, maskVal );
-    
+
     // Smooth the histogram.
     // XXX: is that the right stdev?
     robustHistogramVector = p_psVectorsmoothHistGaussian( robustHistogram,
                             tmpStats->clippedStdev / 4.0f );
-                            
+
     // The following was necessary to fit a gaussian to the data, since
     // gaussian functions produce data between 0.0 and 1.0.
     // p_psNormalizeVector(robustHistogramVector);
-    
+
     /**************************************************************************
     Determine the lower/upper quartiles.
-    **************************************************************************/ 
+    **************************************************************************/
     // We define a vector called "cumulativeRobustSums..." where the value at
     // index position i is equal to the sum of bins 0:i.  This will be used
@@ -1258 +1261 @@
     cumulativeRobustSumsFullRange->data.F32[ 0 ] = robustHistogramVector->data.F32[ 0 ];
     for ( i = 1;i < robustHistogramVector->n;i++ ) {
-            cumulativeRobustSumsFullRange->data.F32[ i ] =
-                cumulativeRobustSumsFullRange->data.F32[ i - 1 ] +
-                robustHistogramVector->data.F32[ i ];
-        }
+        cumulativeRobustSumsFullRange->data.F32[ i ] =
+            cumulativeRobustSumsFullRange->data.F32[ i - 1 ] +
+            robustHistogramVector->data.F32[ i ];
+    }
     sumRobust = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ];
-    
+
     // Determine the bin number containing the lower quartile point.
     LQBinNum = -1;
     for ( i = 0;i < cumulativeRobustSumsFullRange->n;i++ ) {
-            if ( cumulativeRobustSumsFullRange->data.F32[ i ] >= ( sumRobust / 4.0 ) ) {
-                    LQBinNum = i;
-                    break;
-                }
-        }
-        
+        if ( cumulativeRobustSumsFullRange->data.F32[ i ] >= ( sumRobust / 4.0 ) ) {
+            LQBinNum = i;
+            break;
+        }
+    }
+
     // Determine the bin number containing the upper quartile point.
     UQBinNum = -1;
     for ( i = cumulativeRobustSumsFullRange->n - 1;i >= 0;i-- ) {
-            if ( cumulativeRobustSumsFullRange->data.F32[ i ] <= ( 3.0 * sumRobust / 4.0 ) ) {
-                    UQBinNum = i;
-                    break;
-                }
-        }
-        
+        if ( cumulativeRobustSumsFullRange->data.F32[ i ] <= ( 3.0 * sumRobust / 4.0 ) ) {
+            UQBinNum = i;
+            break;
+        }
+    }
+
     if ( ( LQBinNum == -1 ) ||
             ( UQBinNum == -1 ) ) {
-            psAbort( __func__, "Could not determine the robust lower/upper quartiles." );
-        }
+        psAbort( __func__, "Could not determine the robust lower/upper quartiles." );
+    }
     /**************************************************************************
     Determine the mode in the range LQ:UQ.
-    **************************************************************************/ 
+    **************************************************************************/
     // Determine the bin with the peak value in the range LQ to UQ.
     maxBinNum = LQBinNum;
@@ -1294 +1297 @@
     sumN50 = ( float ) robustHistogram->nums->data.S32[ LQBinNum ];
     for ( i = LQBinNum + 1;i <= UQBinNum;i++ ) {
-            if ( robustHistogramVector->data.F32[ i ] > maxBinCount ) {
-                    maxBinNum = i;
-                    maxBinCount = robustHistogramVector->data.F32[ i ];
-                }
-            sumN50 += ( float ) robustHistogram->nums->data.S32[ i ];
-        }
-        
+        if ( robustHistogramVector->data.F32[ i ] > maxBinCount ) {
+            maxBinNum = i;
+            maxBinCount = robustHistogramVector->data.F32[ i ];
+        }
+        sumN50 += ( float ) robustHistogram->nums->data.S32[ i ];
+    }
+
     // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
     dL = ( UQBinNum - LQBinNum ) / 4;
-    
+
     printf( "(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum );
-    
+
     /**************************************************************************
     Determine the mean/stdev for the bins in the range mode-dL to mode+dL
@@ -1311 +1314 @@
     cumulativeRobustSumsDlRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 );
     for ( i = 0;i < robustHistogramVector->n;i++ ) {
-            cumulativeRobustSumsDlRange->data.F32[ i ] = 0.0;
-        }
+        cumulativeRobustSumsDlRange->data.F32[ i ] = 0.0;
+    }
     sumNfit = 0.0;
     cumulativeMedian = 0.0;
     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
-            if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
-                    cumulativeRobustSumsDlRange->data.F32[ i ] =
-                        cumulativeRobustSumsDlRange->data.F32[ i - 1 ] +
-                        robustHistogramVector->data.F32[ i ];
-                    cumulativeMedian += robustHistogramVector->data.F32[ i ];
-                    sumNfit += ( float ) robustHistogram->nums->data.S32[ i ];
-                }
-        }
-        
+        if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
+            cumulativeRobustSumsDlRange->data.F32[ i ] =
+                cumulativeRobustSumsDlRange->data.F32[ i - 1 ] +
+                robustHistogramVector->data.F32[ i ];
+            cumulativeMedian += robustHistogramVector->data.F32[ i ];
+            sumNfit += ( float ) robustHistogram->nums->data.S32[ i ];
+        }
+    }
+
     // Calculate the mean of the smoothed robust histogram in the range
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
@@ -1330 +1333 @@
     myMean = 0.0;
     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
-            if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
-                    myMean += ( robustHistogramVector->data.F32[ i ] ) * 0.5 *
-                              ( robustHistogram->bounds->data.F32[ i + 1 ] +
-                                robustHistogram->bounds->data.F32[ i ] );
-                    countFloat += robustHistogramVector->data.F32[ i ];
-                }
-        }
+        if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
+            myMean += ( robustHistogramVector->data.F32[ i ] ) * 0.5 *
+                      ( robustHistogram->bounds->data.F32[ i + 1 ] +
+                        robustHistogram->bounds->data.F32[ i ] );
+            countFloat += robustHistogramVector->data.F32[ i ];
+        }
+    }
     myMean /= countFloat;
-    
+
     // Calculate the stdev of the smoothed robust histogram in the range
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin.
     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
-            if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
-                    diff = ( 0.5 * ( robustHistogram->bounds->data.F32[ i + 1 ] +
-                                     robustHistogram->bounds->data.F32[ i ] ) ) - myMean;
-                    sumSquares += diff * diff * robustHistogramVector->data.F32[ i ];
-                    sumDiffs += diff * robustHistogramVector->data.F32[ i ];
-                }
-        }
+        if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
+            diff = ( 0.5 * ( robustHistogram->bounds->data.F32[ i + 1 ] +
+                             robustHistogram->bounds->data.F32[ i ] ) ) - myMean;
+            sumSquares += diff * diff * robustHistogramVector->data.F32[ i ];
+            sumDiffs += diff * robustHistogramVector->data.F32[ i ];
+        }
+    }
     myStdev = sqrt( ( sumSquares - ( sumDiffs * sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
-    
+
     /**************************************************************************
     Set the appropriate members in the output stats struct.
     **************************************************************************/
     if ( stats->options & PS_STAT_ROBUST_MEAN ) {
-            stats->robustMean = myMean;
-        }
-        
+        stats->robustMean = myMean;
+    }
+
     if ( stats->options & PS_STAT_ROBUST_MODE ) {
-            stats->robustMode = 0.5 *
-                                ( robustHistogram->bounds->data.F32[ maxBinNum ] +
-                                  robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
-        }
-        
+        stats->robustMode = 0.5 *
+                            ( robustHistogram->bounds->data.F32[ maxBinNum ] +
+                              robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
+    }
+
     if ( stats->options & PS_STAT_ROBUST_STDEV ) {
-            stats->robustStdev = myStdev;
-        }
-        
+        stats->robustStdev = myStdev;
+    }
+
     // To determine the median (and later, the lower/upper quartile), we fit
     // a quadratic to the three bins surrounding the bin containing the median.
@@ -1374 +1377 @@
     // the cumulative number of data points in all bins up to, and including,
     // this bin.  We then solve the quadratic for
-    
+
     if ( stats->options & PS_STAT_ROBUST_MEDIAN ) {
-            if ( ( maxBinNum > 0 ) && ( maxBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
-                    x->data.F64[ 0 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ maxBinNum - 1 ] +
-                                         robustHistogram->bounds->data.F32[ maxBinNum ] );
-                    x->data.F64[ 1 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ maxBinNum ] +
-                                         robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
-                    x->data.F64[ 2 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
-                                         robustHistogram->bounds->data.F32[ maxBinNum + 2 ] );
-                                         
-                    y->data.F64[ 0 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum - 1 ];
-                    y->data.F64[ 1 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum ];
-                    y->data.F64[ 2 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum + 1 ];
-                    
-                    // Ensure that cumulativeMedian/2 is actually within the range of the bins
-                    // we are using.
-                    cumulativeMedian *= 0.5;
-                    if ( !( ( y->data.F64[ 0 ] <= cumulativeMedian ) &&
-                            ( cumulativeMedian <= y->data.F64[ 2 ] ) ) ) {
-                            printf( "((%f), %f, %f)\n", cumulativeMedian, y->data.F64[ 0 ], y->data.F64[ 2 ] );
-                            psAbort( __func__, "p_psVectorRobustStats(1): midpoint not within y-range\n" );
-                        }
-                    // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
-                    // may have to set this meaningfully later.
-                    yErr->data.F64[ 0 ] = 1.0;
-                    yErr->data.F64[ 1 ] = 1.0;
-                    yErr->data.F64[ 2 ] = 1.0;
-                    
-                    // Determine the coefficients of the polynomial.
-                    myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
-                    // Call p_ps1DPolyMedian(), which does a binary search on the
-                    // polynomial, looking for the value x such that
-                    // f(x) = cumulativeMedian.
-                    stats->robustMedian = p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ],
-                                                            x->data.F64[ 2 ], cumulativeMedian );
-                } else {
-                    // If the mode is the first/last histogram bin, then simply use
-                    // the midpoint of that bin.
-                    stats->robustMedian = 0.5 * ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
-                                                  robustHistogram->bounds->data.F32[ maxBinNum ] );
-                }
-        }
-        
+        if ( ( maxBinNum > 0 ) && ( maxBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
+            x->data.F64[ 0 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ maxBinNum - 1 ] +
+                                 robustHistogram->bounds->data.F32[ maxBinNum ] );
+            x->data.F64[ 1 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ maxBinNum ] +
+                                 robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
+            x->data.F64[ 2 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
+                                 robustHistogram->bounds->data.F32[ maxBinNum + 2 ] );
+
+            y->data.F64[ 0 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum - 1 ];
+            y->data.F64[ 1 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum ];
+            y->data.F64[ 2 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum + 1 ];
+
+            // Ensure that cumulativeMedian/2 is actually within the range of the bins
+            // we are using.
+            cumulativeMedian *= 0.5;
+            if ( !( ( y->data.F64[ 0 ] <= cumulativeMedian ) &&
+                    ( cumulativeMedian <= y->data.F64[ 2 ] ) ) ) {
+                printf( "((%f), %f, %f)\n", cumulativeMedian, y->data.F64[ 0 ], y->data.F64[ 2 ] );
+                psAbort( __func__, "p_psVectorRobustStats(1): midpoint not within y-range\n" );
+            }
+            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
+            // may have to set this meaningfully later.
+            yErr->data.F64[ 0 ] = 1.0;
+            yErr->data.F64[ 1 ] = 1.0;
+            yErr->data.F64[ 2 ] = 1.0;
+
+            // Determine the coefficients of the polynomial.
+            myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
+            // Call p_ps1DPolyMedian(), which does a binary search on the
+            // polynomial, looking for the value x such that
+            // f(x) = cumulativeMedian.
+            stats->robustMedian = p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ],
+                                                    x->data.F64[ 2 ], cumulativeMedian );
+        } else {
+            // If the mode is the first/last histogram bin, then simply use
+            // the midpoint of that bin.
+            stats->robustMedian = 0.5 * ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
+                                          robustHistogram->bounds->data.F32[ maxBinNum ] );
+        }
+    }
+
     // The lower/upper quartile calculations are very similar to the median
     // calculations.  We fit a quadratic to the array containing the
@@ -1426 +1429 @@
     //
     if ( stats->options & PS_STAT_ROBUST_QUARTILE ) {
-            countFloat = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ];
-            
-            if ( ( LQBinNum > 0 ) && ( LQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
-                    x->data.F64[ 0 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ LQBinNum - 1 ] +
-                                         robustHistogram->bounds->data.F32[ LQBinNum ] );
-                    x->data.F64[ 1 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ LQBinNum ] +
-                                         robustHistogram->bounds->data.F32[ LQBinNum + 1 ] );
-                    x->data.F64[ 2 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
-                                         robustHistogram->bounds->data.F32[ LQBinNum + 2 ] );
-                                         
-                    y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum - 1 ];
-                    y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum ];
-                    y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum + 1 ];
-                    
-                    if ( !( ( y->data.F64[ 0 ] <= ( countFloat / 4.0 ) ) &&
-                            ( ( countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
-                            psAbort( __func__, "p_psVectorRobustStats(2): midpoint not within y-range\n" );
-                        }
-                        
-                    yErr->data.F64[ 0 ] = 1.0;
-                    yErr->data.F64[ 1 ] = 1.0;
-                    yErr->data.F64[ 2 ] = 1.0;
-                    
-                    myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
-                    stats->robustLQ = p_ps1DPolyMedian( myPoly,
-                                                        x->data.F64[ 0 ],
-                                                        x->data.F64[ 2 ],
-                                                        countFloat / 4.0 );
-                                                        
-                } else {
-                    // If the LQ is the first/last histogram bin, then simply use
-                    // the midpoint of that bin.
-                    stats->robustLQ = 0.5 * ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
-                                              robustHistogram->bounds->data.F32[ LQBinNum ] );
-                }
-                
-            if ( ( UQBinNum > 0 ) && ( UQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
-                    x->data.F64[ 0 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ UQBinNum - 1 ] +
-                                         robustHistogram->bounds->data.F32[ UQBinNum ] );
-                    x->data.F64[ 1 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ UQBinNum ] +
-                                         robustHistogram->bounds->data.F32[ UQBinNum + 1 ] );
-                    x->data.F64[ 2 ] = ( double ) 0.5 *
-                                       ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
-                                         robustHistogram->bounds->data.F32[ UQBinNum + 2 ] );
-                                         
-                    y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum - 1 ];
-                    y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum ];
-                    y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum + 1 ];
-                    
-                    if ( !( ( y->data.F64[ 0 ] <= ( 3.0 * countFloat / 4.0 ) ) &&
-                            ( ( 3.0 * countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
-                            psAbort( __func__, "p_psVectorRobustStats(3): midpoint not within y-range\n" );
-                        }
-                        
-                    yErr->data.F64[ 0 ] = 1.0;
-                    yErr->data.F64[ 1 ] = 1.0;
-                    yErr->data.F64[ 2 ] = 1.0;
-                    
-                    myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
-                    stats->robustUQ = p_ps1DPolyMedian( myPoly,
-                                                        x->data.F64[ 0 ],
-                                                        x->data.F64[ 2 ],
-                                                        3.0 * countFloat / 4.0 );
-                } else {
-                    // If the UQ is the first/last histogram bin, then simply use
-                    // the midpoint of that bin.
-                    stats->robustUQ = 0.5 * ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
-                                              robustHistogram->bounds->data.F32[ UQBinNum ] );
-                }
-        }
+        countFloat = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ];
+
+        if ( ( LQBinNum > 0 ) && ( LQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
+            x->data.F64[ 0 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ LQBinNum - 1 ] +
+                                 robustHistogram->bounds->data.F32[ LQBinNum ] );
+            x->data.F64[ 1 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ LQBinNum ] +
+                                 robustHistogram->bounds->data.F32[ LQBinNum + 1 ] );
+            x->data.F64[ 2 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
+                                 robustHistogram->bounds->data.F32[ LQBinNum + 2 ] );
+
+            y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum - 1 ];
+            y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum ];
+            y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum + 1 ];
+
+            if ( !( ( y->data.F64[ 0 ] <= ( countFloat / 4.0 ) ) &&
+                    ( ( countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
+                psAbort( __func__, "p_psVectorRobustStats(2): midpoint not within y-range\n" );
+            }
+
+            yErr->data.F64[ 0 ] = 1.0;
+            yErr->data.F64[ 1 ] = 1.0;
+            yErr->data.F64[ 2 ] = 1.0;
+
+            myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
+            stats->robustLQ = p_ps1DPolyMedian( myPoly,
+                                                x->data.F64[ 0 ],
+                                                x->data.F64[ 2 ],
+                                                countFloat / 4.0 );
+
+        } else {
+            // If the LQ is the first/last histogram bin, then simply use
+            // the midpoint of that bin.
+            stats->robustLQ = 0.5 * ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
+                                      robustHistogram->bounds->data.F32[ LQBinNum ] );
+        }
+
+        if ( ( UQBinNum > 0 ) && ( UQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
+            x->data.F64[ 0 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ UQBinNum - 1 ] +
+                                 robustHistogram->bounds->data.F32[ UQBinNum ] );
+            x->data.F64[ 1 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ UQBinNum ] +
+                                 robustHistogram->bounds->data.F32[ UQBinNum + 1 ] );
+            x->data.F64[ 2 ] = ( double ) 0.5 *
+                               ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
+                                 robustHistogram->bounds->data.F32[ UQBinNum + 2 ] );
+
+            y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum - 1 ];
+            y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum ];
+            y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum + 1 ];
+
+            if ( !( ( y->data.F64[ 0 ] <= ( 3.0 * countFloat / 4.0 ) ) &&
+                    ( ( 3.0 * countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
+                psAbort( __func__, "p_psVectorRobustStats(3): midpoint not within y-range\n" );
+            }
+
+            yErr->data.F64[ 0 ] = 1.0;
+            yErr->data.F64[ 1 ] = 1.0;
+            yErr->data.F64[ 2 ] = 1.0;
+
+            myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
+            stats->robustUQ = p_ps1DPolyMedian( myPoly,
+                                                x->data.F64[ 0 ],
+                                                x->data.F64[ 2 ],
+                                                3.0 * countFloat / 4.0 );
+        } else {
+            // If the UQ is the first/last histogram bin, then simply use
+            // the midpoint of that bin.
+            stats->robustUQ = 0.5 * ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
+                                      robustHistogram->bounds->data.F32[ UQBinNum ] );
+        }
+    }
     stats->robustNfit = sumNfit;
     stats->robustN50 = sumN50;
-    
+
     psFree( x );
     psFree( y );
@@ -1576 +1579 @@
 {
     psStats * newStruct = NULL;
-    
+
     newStruct = ( psStats * ) psAlloc( sizeof( psStats ) );
     newStruct->sampleMean = NAN;
@@ -1600 +1603 @@
     newStruct->binsize = NAN;
     newStruct->options = options;
-    
+
     return ( newStruct );
 }
@@ -1623 +1626 @@
     psHistogram *newHist = NULL;        // The new histogram structure
     float binSize = 0.0;                // The histogram bin size
-    
+
     // NOTE: Verify that this is the correct action.
     if ( n == 0 ) {
-            return ( NULL );
-        }
+        return ( NULL );
+    }
     if ( n < 0 ) {
-            psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n );
-        }
-        
+        psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n );
+    }
+
     // NOTE: Verify that this is the correct action.
     if ( lower > upper ) {
-            return ( NULL );
-        }
-        
+        return ( NULL );
+    }
+
     // Allocate memory for the new histogram structure.  If there are N
     // bins, then there are N+1 bounds to those bins.
@@ -1643 +1646 @@
     newHist->bounds = psVectorAlloc( n + 1, PS_TYPE_F32 );
     newHist->bounds->n = newHist->bounds->nalloc;
-    
+
     // Calculate the bounds for each bin.
     binSize = ( upper - lower ) / ( float ) n;
@@ -1650 +1653 @@
     binSize += FLT_EPSILON;
     for ( i = 0;i < n + 1;i++ ) {
-            newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) i );
-        }
-        
+        newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) i );
+    }
+
     // Allocate the bins, and initialize them to zero.
     newHist->nums = psVectorAlloc( n, PS_TYPE_U32 );
     newHist->nums->n = newHist->nums->nalloc;
     for ( i = 0;i < newHist->nums->n;i++ ) {
-            newHist->nums->data.U32[ i ] = 0;
-        }
-        
+        newHist->nums->data.U32[ i ] = 0;
+    }
+
     // Initialize the other members.
     newHist->minNum = 0;
     newHist->maxNum = 0;
     newHist->uniform = true;
-    
+
     return ( newHist );
 }
@@ -1681 +1684 @@
     psHistogram * newHist = NULL;        // The new histogram structure
     int i;                              // Loop index variable
-    
+
     // NOTE: Verify that this is the correct action.
     if ( bounds == NULL ) {
-            // psAbort(__func__, "psHistogram requested with NULL bounds");
-            return ( NULL );
-        }
-        
+        // psAbort(__func__, "psHistogram requested with NULL bounds");
+        return ( NULL );
+    }
+
     // NOTE: Verify that this is the correct action.
     if ( bounds->n <= 1 ) {
-            // psAbort(__func__, "psHistogram requested with NULL bounds");
-            return ( NULL );
-        }
-        
+        // psAbort(__func__, "psHistogram requested with NULL bounds");
+        return ( NULL );
+    }
+
     if ( bounds->type.type != PS_TYPE_F32 ) {
-            // psAbort(__func__, "psHistogram request a bound which is not type F32");
-            return ( NULL );
-        }
-        
+        // psAbort(__func__, "psHistogram request a bound which is not type F32");
+        return ( NULL );
+    }
+
     // Allocate memory for the new histogram structure.
     newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) );
@@ -1705 +1708 @@
     newHist->bounds->n = newHist->bounds->nalloc;
     for ( i = 0;i < bounds->n;i++ ) {
-            newHist->bounds->data.F32[ i ] = bounds->data.F32[ i ];
-        }
-        
+        newHist->bounds->data.F32[ i ] = bounds->data.F32[ i ];
+    }
+
     // Allocate the bins, and initialize them to zero.  If there are N bounds,
     // then there are N-1 bins.
@@ -1713 +1716 @@
     newHist->nums->n = newHist->nums->nalloc;
     for ( i = 0;i < newHist->nums->n;i++ ) {
-            newHist->nums->data.U32[ i ] = 0;
-        }
-        
+        newHist->nums->data.U32[ i ] = 0;
+    }
+
     // Initialize the other members.
     newHist->minNum = 0;
     newHist->maxNum = 0;
     newHist->uniform = false;
-    
+
     return ( newHist );
 }
@@ -1755 +1758 @@
     int binNum = 0;                             // A temporary bin number
     int numBins = 0;                            // The total number of bins
-    
+
     // NOTE: Verify that this is the correct action.
     if ( out == NULL ) {
-            return ( NULL );
-        }
-        
+        return ( NULL );
+    }
+
     // Check the specified output histogram for type psF32
     if ( out->bounds->type.type != PS_TYPE_F32 ) {
-            psAbort( __func__,
-                     "Only data type PS_TYPE_F32 for the output.bounds member." );
-        }
-        
+        psAbort( __func__,
+                 "Only data type PS_TYPE_F32 for the output.bounds member." );
+    }
+
     if ( out->nums->type.type != PS_TYPE_U32 ) {
-            psAbort( __func__,
-                     "Only data type PS_TYPE_U32 for output.nums member." );
-        }
-        
+        psAbort( __func__,
+                 "Only data type PS_TYPE_U32 for output.nums member." );
+    }
+
     // NOTE: Verify that this is the correct action.
     if ( in == NULL ) {
-            return ( out );
-        }
-        
+        return ( out );
+    }
+
     if ( in->type.type != PS_TYPE_F32 ) {
+        psAbort( __func__,
+                 "Only data type PS_TYPE_F32 is currently supported (0x%x).",
+                 in->type.type );
+    }
+
+    if ( mask != NULL ) {
+        if ( in->n != mask->n ) {
             psAbort( __func__,
-                     "Only data type PS_TYPE_F32 is currently supported (0x%x).",
-                     in->type.type );
-        }
-        
-    if ( mask != NULL ) {
-            if ( in->n != mask->n ) {
-                    psAbort( __func__,
-                             "Vector data and vector mask are of different sizes." );
-                }
-            if ( mask->type.type != PS_TYPE_U8 ) {
-                    psAbort( __func__, "Vector mask must be type PS_TYPE_U8" );
-                }
-        }
+                     "Vector data and vector mask are of different sizes." );
+        }
+        if ( mask->type.type != PS_TYPE_U8 ) {
+            psAbort( __func__, "Vector mask must be type PS_TYPE_U8" );
+        }
+    }
     // NOTE: determine the correct action for a variety of other cases:
     // in vector has 0 elements, and histogram structure has zero bins.
-    
+
     numBins = out->nums->n;
     for ( i = 0;i < in->n;i++ ) {
-            // Check if this pixel is masked, and if so, skip it.
-            if ( ( mask == NULL ) ||
-                    ( ( mask != NULL ) && ( !( 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 );
-                                                       
-                                    // NOTE: This next if-statement really shouldn't be necessary.
-                                    // However, do to numerical lack of precision, we occasionally
-                                    // produce a binNum outside the range of bins.
-                                    if ( binNum >= out->nums->n ) {
-                                            binNum = out->nums->n - 1;
-                                        }
-                                        
-                                    ( out->nums->data.S32[ binNum ] ) ++;
-                                    
-                                    // If this is a non-uniform histogram, determining the correct
-                                    // bin number requires a bit more work.
-                                } else {
-                                    // NOTE: 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 ] ) ++;
-                                                }
-                                        }
-                                }
+        // Check if this pixel is masked, and if so, skip it.
+        if ( ( mask == NULL ) ||
+                ( ( mask != NULL ) && ( !( 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 );
+
+                        // NOTE: This next if-statement really shouldn't be necessary.
+                        // However, do to numerical lack of precision, we occasionally
+                        // produce a binNum outside the range of bins.
+                        if ( binNum >= out->nums->n ) {
+                            binNum = out->nums->n - 1;
                         }
-                }
-        }
+
+                        ( out->nums->data.S32[ binNum ] ) ++;
+
+                        // If this is a non-uniform histogram, determining the correct
+                        // bin number requires a bit more work.
+                    } else {
+                        // NOTE: 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 ] ) ++;
+                            }
+                        }
+                    }
+                }
+        }
+    }
     return ( out );
 }
@@ -1859 +1863 @@
     int i = 0;
     psVector *tmp = NULL;
-    
+
     if ( in->type.type == PS_TYPE_S32 ) {
+        tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
+        for ( i = 0;i < in->n;i++ ) {
+            tmp->data.F32[ i ] = ( float ) in->data.S32[ i ];
+        }
+    } else
+        if ( in->type.type == PS_TYPE_U32 ) {
             tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
             for ( i = 0;i < in->n;i++ ) {
-                    tmp->data.F32[ i ] = ( float ) in->data.S32[ i ];
-                }
-        } else if ( in->type.type == PS_TYPE_U32 ) {
-            tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
-            for ( i = 0;i < in->n;i++ ) {
-                    tmp->data.F32[ i ] = ( float ) in->data.U32[ i ];
-                }
-        } else if ( in->type.type == PS_TYPE_F64 ) {
-            tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
-            for ( i = 0;i < in->n;i++ ) {
+                tmp->data.F32[ i ] = ( float ) in->data.U32[ i ];
+            }
+        } else
+            if ( in->type.type == PS_TYPE_F64 ) {
+                tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
+                for ( i = 0;i < in->n;i++ ) {
                     tmp->data.F32[ i ] = ( float ) in->data.F64[ i ];
                 }
-        } else if ( in->type.type == PS_TYPE_F32 ) {
-            // do nothing
-        } else {
-            psAbort( __func__, "unsupported vector type 0x%x\n", in->type.type );
-        }
+            } else
+                if ( in->type.type == PS_TYPE_F32 ) {
+                    // do nothing
+                } else {
+                    psAbort( __func__, "unsupported vector type 0x%x\n", in->type.type );
+                }
     return ( tmp );
 }
@@ -1908 +1915 @@
     psVector * inF32;
     int mustFreeTmp = 1;
-    
+
     // NOTE: Verify that this is the correct action.
     if ( in == NULL ) {
-            return ( stats );
-        }
+        return ( stats );
+    }
     if ( stats == NULL ) {
-            return ( NULL );
-        }
-        
+        return ( NULL );
+    }
+
     inF32 = p_psConvertToF32( stats, in, mask, maskVal );
     if ( inF32 == NULL ) {
-            inF32 = in;
-            mustFreeTmp = 0;
-        }
-        
+        inF32 = in;
+        mustFreeTmp = 0;
+    }
+
     // XXX: Should we abort if (stats->min == stats->max)?
     if ( ( stats->options & PS_STAT_USE_RANGE ) &&
             ( stats->min >= stats->max ) ) {
-            psAbort( __func__, "psVectorStats() called with range: %f to %f\n",
-                     stats->min, stats->max );
-        }
-        
+        psAbort( __func__, "psVectorStats() called with range: %f to %f\n",
+                 stats->min, stats->max );
+    }
+
     //    PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
     if ( mask != NULL ) {
-            PS_CHECK_NULL_VECTOR( mask );
-            PS_CHECK_EMPTY_VECTOR( mask );
-            PS_CHECK_VECTOR_SIZE_EQUAL( mask, in );
-            PS_CHECK_VECTOR_TYPE( mask, PS_TYPE_U8 );
-        }
-        
+        PS_CHECK_NULL_VECTOR( mask );
+        PS_CHECK_EMPTY_VECTOR( mask );
+        PS_CHECK_VECTOR_SIZE_EQUAL( mask, in );
+        PS_CHECK_VECTOR_TYPE( mask, PS_TYPE_U8 );
+    }
+
     // ************************************************************************
     if ( stats->options & PS_STAT_SAMPLE_MEAN ) {
-            p_psVectorSampleMean( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorSampleMean( in, mask, maskVal, stats );
+    }
+
     // ************************************************************************
     if ( stats->options & PS_STAT_SAMPLE_MEDIAN ) {
-            p_psVectorSampleMedian( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorSampleMedian( in, mask, maskVal, stats );
+    }
+
     // ************************************************************************
     // NOTE: 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 );
-        }
-        
+        p_psVectorSampleMean( in, mask, maskVal, stats );
+        p_psVectorSampleStdev( in, mask, maskVal, stats );
+    }
+
     // ************************************************************************
     if ( stats->options & PS_STAT_SAMPLE_QUARTILE ) {
-            p_psVectorSampleQuartiles( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorSampleQuartiles( in, mask, maskVal, stats );
+    }
+
     // Since the various robust stats quantities share much computation, they
     // are grouped together in a single private function:
@@ -1969 +1976 @@
             ( stats->options & PS_STAT_ROBUST_STDEV ) ||
             ( stats->options & PS_STAT_ROBUST_QUARTILE ) ) {
-            p_psVectorRobustStats( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorRobustStats( in, mask, maskVal, stats );
+    }
+
     if ( ( stats->options & PS_STAT_CLIPPED_MEAN ) ||
             ( stats->options & PS_STAT_CLIPPED_STDEV ) ) {
-            p_psVectorClippedStats( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorClippedStats( in, mask, maskVal, stats );
+    }
+
     // ************************************************************************
     if ( stats->options & PS_STAT_MAX ) {
-            p_psVectorMax( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorMax( in, mask, maskVal, stats );
+    }
+
     // ************************************************************************
     if ( stats->options & PS_STAT_MIN ) {
-            p_psVectorMin( in, mask, maskVal, stats );
-        }
-        
+        p_psVectorMin( in, mask, maskVal, stats );
+    }
+
     if ( mustFreeTmp == 1 ) {
-            psFree( inF32 );
-        }
+        psFree( inF32 );
+    }
     return ( stats );
 }