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

Timestamp:

Jul 31, 2004, 3:16:01 PM (22 years ago)

Author:

gusciora

Message:

...

File:

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

trunk/psLib/src/math/psStats.c

@@ -32 +32 @@
 #define MAX_ITERATIONS 10
 
-void p_psVectorRobustStats(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *stats);
-
+void p_psVectorRobustStats( const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            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"); \
 
@@ -87 +87 @@
 /*****************************************************************************/
 
-bool p_psGetStatValue(const psStats* stats, double* value)
-{
-
-    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;
-    }
+bool p_psGetStatValue( const psStats* stats, double* value )
+{
+
+    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;
+        }
 }
 
@@ -147 +147 @@
  *****************************************************************************/
 
-void p_psVectorPrint(psVector *myVector,
-                     psVector *maskVector,
-                     unsigned int maskVal,
-                     psStats *stats)
+void p_psVectorPrint( psVector *myVector,
+                      psVector *maskVector,
+                      unsigned int maskVal,
+                      psStats *stats )
 {
     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]);
-    }
+    
+    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 ] );
+        }
 }
 
@@ -194 +194 @@
  *****************************************************************************/
 
-void p_psVectorSampleMean(const psVector *restrict myVector,
-                          const psVector *restrict maskVector,
-                          unsigned int maskVal,
-                          psStats *stats)
+void p_psVectorSampleMean( const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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;
+    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;
+                }
         } 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;
-        }
-    } 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;
-        }
-    }
-
+            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;
 }
@@ -262 +262 @@
     NULL
  *****************************************************************************/
-void p_psVectorMax(const psVector *restrict myVector,
-                   const psVector *restrict maskVector,
-                   unsigned int maskVal,
-                   psStats *stats)
+void p_psVectorMax( const psVector *restrict myVector,
+                    const psVector *restrict maskVector,
+                    unsigned int maskVal,
+                    psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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];
-                    }
-                }
-            }
+    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 ];
+                                }
+                        }
+                }
         } 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];
-                }
-            }
-        }
-    } 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];
-                }
-            }
-        }
-    }
-
+            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;
 }
@@ -327 +327 @@
     NULL
  *****************************************************************************/
-void p_psVectorMin(const psVector *restrict myVector,
-                   const psVector *restrict maskVector,
-                   unsigned int maskVal,
-                   psStats *stats)
+void p_psVectorMin( const psVector *restrict myVector,
+                    const psVector *restrict maskVector,
+                    unsigned int maskVal,
+                    psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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];
-                    }
-                }
-            }
+    
+    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 ];
+                                }
+                        }
+                }
         } 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];
-                }
-            }
-        }
-    } 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];
-                }
-            }
-        }
-    }
-
+            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;
 }
@@ -392 +392 @@
     NULL
  *****************************************************************************/
-int p_psVectorNValues(const psVector *restrict myVector,
-                      const psVector *restrict maskVector,
-                      unsigned int maskVal,
-                      psStats *stats)
+int p_psVectorNValues( const psVector *restrict myVector,
+                       const psVector *restrict maskVector,
+                       unsigned int maskVal,
+                       psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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++;
-                }
-            }
+    
+    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++;
+                                }
+                        }
+                }
         } else {
-            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);
+            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 );
 }
 
@@ -450 +450 @@
     NULL
  *****************************************************************************/
-void p_psVectorSampleMedian(const psVector *restrict myVector,
-                            const psVector *restrict maskVector,
-                            unsigned int maskVal,
-                            psStats *stats)
-{
-    psVector *unsortedVector = NULL;            // Temporary vector
+void p_psVectorSampleMedian( const psVector *restrict myVector,
+                             const psVector *restrict maskVector,
+                             unsigned int maskVal,
+                             psStats *stats )
+{
+    psVector * unsortedVector = NULL;            // Temporary vector
     psVector *sortedVector = NULL;              // Temporary vector
     int i = 0;                                  // Loop index variable
@@ -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);
-
+    nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats );
+    
     // Allocate temporary vectors for the data.
-    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
+    unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     unsortedVector->n = unsortedVector->nalloc;
-
-    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
+    
+    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];
-                }
-            }
+    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 ];
+                                }
+                        }
+                }
         } 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];
-                }
-            }
-        }
-    } 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];
-            }
-        }
-    }
+            // 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);
-
+    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];
-    }
-
+    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 ];
+        }
+        
     // Free the temporary data structures.
-    psFree(unsortedVector);
-    psFree(sortedVector);
+    psFree( unsortedVector );
+    psFree( sortedVector );
 }
 
@@ -568 +568 @@
     XXX: use a static variable for gaussianCoefs[] and compute them once.
  *****************************************************************************/
-psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
-                                       float sigma)
+psVector *p_psVectorsmoothHistGaussian( psHistogram *robustHistogram,
+                                        float sigma )
 {
     int i = 0;                                  // Loop index variable
@@ -575 +575 @@
     float denom = 0.0;                           // Temporary variable
     float expo = 0.0;                            // Temporary variable
-    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);
-        }
-    }
-
+    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 );
+                }
+        }
+        
     // 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]);
-            }
-        }
-    }
-    return(smooth);
+    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 ] );
+                        }
+                }
+        }
+    return ( smooth );
 }
 
@@ -634 +634 @@
     NULL
  *****************************************************************************/
-void p_psVectorSampleQuartiles(const psVector *restrict myVector,
-                               const psVector *restrict maskVector,
-                               unsigned int maskVal,
-                               psStats *stats)
-{
-    psVector *unsortedVector = NULL;            // Temporary vector
+void p_psVectorSampleQuartiles( const psVector *restrict myVector,
+                                const psVector *restrict maskVector,
+                                unsigned int maskVal,
+                                psStats *stats )
+{
+    psVector * unsortedVector = NULL;            // Temporary vector
     psVector *sortedVector = NULL;              // Temporary vector
     int i = 0;                                  // Loop index variable
@@ -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);
-
+    nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats );
+    
     // Allocate temporary vectors for the data.
-    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
+    unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     unsortedVector->n = unsortedVector->nalloc;
-    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
+    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];
-                }
-            }
+    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 ];
+                                }
+                        }
+                }
         } 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];
-                }
-            }
-        }
-    } 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];
-            }
-        }
-    }
-
+            // 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);
-
+    psVectorSort( sortedVector, unsortedVector );
+    
     // Calculate the quartile points exactly.
-    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
-    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
-
+    stats->sampleUQ = sortedVector->data.F32[ 3 * ( nValues / 4 ) ];
+    stats->sampleLQ = sortedVector->data.F32[ nValues / 4 ];
+    
     // Free the temporary data structures.
-    psFree(unsortedVector);
-    psFree(sortedVector);
+    psFree( unsortedVector );
+    psFree( sortedVector );
 }
 
@@ -721 +721 @@
  
  *****************************************************************************/
-void p_psVectorSampleStdev(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *stats)
+void p_psVectorSampleStdev( const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal,
+                            psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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);
-    }
+    if ( 0 != isnan( stats->sampleMean ) ) {
+            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++;
-                }
-            }
+    
+    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;
+                }
         } 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;
-        }
-    } 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;
-
+            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));
+    
+    stats->sampleStdev = ( float ) sqrt( ( sumSquares - ( sumDiffs *
+                                           sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     #else
-
-    stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
-                                 sumDiffs/countFloat))/ (countFloat-1));
+    
+    stats->sampleStdev = sqrtf( ( sumSquares - ( sumDiffs *
+                                  sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     #endif
 }
@@ -812 +812 @@
     NULL
  *****************************************************************************/
-void p_psVectorClippedStats(const psVector *restrict myVector,
-                            const psVector *restrict maskVector,
-                            unsigned int maskVal,
-                            psStats *stats)
+void p_psVectorClippedStats( const psVector *restrict myVector,
+                             const psVector *restrict maskVector,
+                             unsigned int maskVal,
+                             psStats *stats )
 {
     int i = 0;                                  // Loop index variable
@@ -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);
-    }
-
+    if ( !( ( CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
+            ( stats->clipIter <= CLIPPED_NUM_ITER_UB ) ) ) {
+            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);
-    }
-
+    if ( !( ( CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
+            ( stats->clipSigma <= CLIPPED_SIGMA_UB ) ) ) {
+            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.
     // However, we do no want to modify the original mask vector.
-    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+    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];
-        }
-    }
-
+    if ( maskVector != NULL ) {
+            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);
-
+    p_psVectorSampleMedian( myVector, maskVector, maskVal, stats );
+    
     // 2. Compute the sample standard deviation.
-    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
-
+    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 ( 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;
+                }
+                
+        }
     stats->sampleMean = oldStanMean;
-    stats->sampleStdev= oldStanStdev;
-
+    stats->sampleStdev = oldStanStdev;
+    
     // 7. The last calcuated value of x is the cliped mean.
-    if (stats->options & PS_STAT_CLIPPED_MEAN) {
-        stats->clippedMean = clippedMean;
-    }
-
+    if ( stats->options & PS_STAT_CLIPPED_MEAN ) {
+            stats->clippedMean = clippedMean;
+        }
+        
     // 8. The last calcuated value of stdev is the cliped stdev.
-    if (stats->options & PS_STAT_CLIPPED_STDEV) {
-        stats->clippedStdev = clippedStdev;
-    }
-
-    psFree(tmpMask);
+    if ( stats->options & PS_STAT_CLIPPED_STDEV ) {
+            stats->clippedStdev = clippedStdev;
+        }
+        
+    psFree( tmpMask );
 }
 
@@ -912 +912 @@
 elements of a vector to a range between 0.0 and 1.0.
  *****************************************************************************/
-void p_psNormalizeVector(psVector *myData)
-{
-    float min = (float) HUGE;
-    float max = (float) -HUGE;
+void p_psNormalizeVector( psVector *myData )
+{
+    float min = ( float ) HUGE;
+    float max = ( float ) - HUGE;
     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];
-        }
-    }
-
+    
+    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 ];
+                }
+        }
+        
     range = max - min;
-    for (i=0;i<myData->n;i++) {
-        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
-    }
+    for ( i = 0;i < myData->n;i++ ) {
+            myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range;
+        }
 }
 
@@ -940 +940 @@
 specified data point.
  *****************************************************************************/
-float p_psGaussian(const psVector *restrict myData,
-                   const psVector *restrict myParams)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
-    float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev));
-    tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
-
+float p_psGaussian( const psVector *restrict myData,
+                    const psVector *restrict myParams )
+{
+    float x = myData->data.F32[ 0 ];
+    float mean = myParams->data.F32[ 0 ];
+    float stdev = myParams->data.F32[ 1 ];
+    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);
+    return ( tmp );
 }
 
@@ -957 +957 @@
 calculates the specified partial derivative of the above Gaussian function.
  *****************************************************************************/
-float p_psGaussianDeriv(const psVector *restrict myData,
-                        const psVector *restrict myParams,
-                        int whichParam)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
+float p_psGaussianDeriv( const psVector *restrict myData,
+                         const psVector *restrict myParams,
+                         int whichParam )
+{
+    float x = myData->data.F32[ 0 ];
+    float mean = myParams->data.F32[ 0 ];
+    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 stdev.
-        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
-        tmp/= (stdev * stdev * stdev);
-    }
-    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
-
-    return(tmp);
+    
+    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 stdev.
+            tmp = ( x - mean ) * ( x - mean ) * p_psGaussian( myData, myParams );
+            tmp /= ( stdev * stdev * stdev );
+        }
+    printf( "p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp );
+    
+    return ( tmp );
 }
 
@@ -986 +986 @@
 specified data point.
  *****************************************************************************/
-float p_psQuadratic(const psVector *restrict myParams,
-                    const psVector *restrict myCoords)
-{
-    float x = myCoords->data.F32[0];
-    float A = myParams->data.F32[0];
-    float B = myParams->data.F32[1];
-    float C = myParams->data.F32[2];
+float p_psQuadratic( const psVector *restrict myParams,
+                     const psVector *restrict myCoords )
+{
+    float x = myCoords->data.F32[ 0 ];
+    float A = myParams->data.F32[ 0 ];
+    float B = myParams->data.F32[ 1 ];
+    float C = myParams->data.F32[ 2 ];
     float tmp = 0.0;
-
-    tmp = (A * x * x) + (B * x) + C;
-    return(tmp);
+    
+    tmp = ( A * x * x ) + ( B * x ) + C;
+    return ( tmp );
 }
 
@@ -1003 +1003 @@
 calculates the specified partial derivative of the above quadratic function.
  *****************************************************************************/
-float p_psQuadraticDeriv(const psVector *restrict myParams,
-                         const psVector *restrict myCoords,
-                         int whichParamDeriv)
-{
-    float x = myCoords->data.F32[0];
+float p_psQuadraticDeriv( const psVector *restrict myParams,
+                          const psVector *restrict myCoords,
+                          int whichParamDeriv )
+{
+    float x = myCoords->data.F32[ 0 ];
     float tmp = 0.0;
-
-    if (whichParamDeriv == 0) {
-        tmp = x * x;
-    } else if (whichParamDeriv == 1) {
-        tmp = x;
-    } else if (whichParamDeriv == 2) {
-        tmp = 1.0;
-    }
-
-    return(tmp);
+    
+    if ( whichParamDeriv == 0 ) {
+            tmp = x * x;
+        } else if ( whichParamDeriv == 1 ) {
+            tmp = x;
+        } else if ( whichParamDeriv == 2 ) {
+            tmp = 1.0;
+        }
+        
+    return ( tmp );
 }
 
@@ -1030 +1030 @@
 decreasing within that range.
  *****************************************************************************/
-float p_ps1DPolyMedian(psPolynomial1D *myPoly,
-                       float rangeLow,
-                       float rangeHigh,
-                       float getThisValue)
-{
-    int numIterations=0;
+float p_ps1DPolyMedian( psPolynomial1D *myPoly,
+                        float rangeLow,
+                        float rangeHigh,
+                        float getThisValue )
+{
+    int numIterations = 0;
     float midpoint = 0.0;
     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++;
-    }
-    return(midpoint);
+    
+    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++;
+        }
+    return ( midpoint );
 }
 
@@ -1073 +1073 @@
 XXX: This function is currently not being used.
  *****************************************************************************/
-float p_psFitQuadratic(psHistogram *histogram,
-                       psVector *cumulativeSums,
-                       int binNum,
-                       float fitFloat)
-{
-    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
-    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 = psGetArrayPolynomial(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);
-    psFree(yErr);
-    psFree(myPoly);
-    return(0.0);
+float p_psFitQuadratic( psHistogram *histogram,
+                        psVector *cumulativeSums,
+                        int binNum,
+                        float fitFloat )
+{
+    psVector * x = psVectorAlloc( 3, PS_TYPE_F64 );
+    psVector *y = psVectorAlloc( 3, PS_TYPE_F64 );
+    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 ] ) );
+        }
+        
+    psFree( x );
+    psFree( y );
+    psFree( yErr );
+    psFree( myPoly );
+    return ( 0.0 );
 }
 
@@ -1147 +1147 @@
     NULL
 *****************************************************************************/
-void p_psVectorRobustStats(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *stats)
-{
-    psHistogram *robustHistogram = NULL;
+void p_psVectorRobustStats( const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal,
+                            psStats *stats )
+{
+    psHistogram * robustHistogram = NULL;
     psVector *robustHistogramVector = NULL;
     float binSize = 0.0;                        // Size of the histogram bins
@@ -1162 +1162 @@
     float dL = 0.0;
     int numBins = 0;
-    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
+    psStats *tmpStats = psStatsAlloc( PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN );
     //    psImage *domain;
     //    psVector *errors;
@@ -1177 +1177 @@
     float sumSquares = 0.0;
     float sumDiffs = 0.0;
-    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
+    psVector *x = psVectorAlloc( 3, PS_TYPE_F64 );
+    psVector *y = psVectorAlloc( 3, PS_TYPE_F64 );
+    psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 );
+    psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 );
     psVector *cumulativeRobustSumsFullRange = NULL;
     psVector *cumulativeRobustSumsDlRange = NULL;
@@ -1187 +1187 @@
     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
     // that by 10.0;
-    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
+    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 ( 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 ;
+        }
+        
     // Determine minimum and maximum values in the data vector.
-    if (isnan(stats->min)) {
-        p_psVectorMin(myVector, maskVector, maskVal, stats);
-    }
-    if (isnan(stats->max)) {
-        p_psVectorMax(myVector, maskVector, maskVal, stats);
-    }
-
+    if ( isnan( stats->min ) ) {
+            p_psVectorMin( myVector, maskVector, maskVal, stats );
+        }
+    if ( isnan( stats->max ) ) {
+            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
     // bins, not the binSize.  Also, if we get here, we know that
     // binSize != 0.0.
-    numBins = (int) ((stats->max - stats->min) / binSize);
-    robustHistogram = psHistogramAlloc(stats->min,
-                                       stats->max,
-                                       numBins);
-
+    numBins = ( int ) ( ( stats->max - stats->min ) / binSize );
+    robustHistogram = psHistogramAlloc( stats->min,
+                                        stats->max,
+                                        numBins );
+                                        
     // Populate the histogram array.
-    psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
-
+    psVectorHistogram( robustHistogram, myVector, maskVector, maskVal );
+    
     // Smooth the histogram.
     // XXX: is that the right stdev?
-    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
-                            tmpStats->clippedStdev/4.0f);
-
+    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
     // now and later in determining the lower/upper quartiles.
-    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    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];
-    }
-    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1];
-
+    cumulativeRobustSumsFullRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 );
+    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 ];
+        }
+    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;
-        }
-    }
-
+    for ( i = 0;i < cumulativeRobustSumsFullRange->n;i++ ) {
+            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 ((LQBinNum == -1) ||
-            (UQBinNum == -1)) {
-        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
-    }
+    for ( i = cumulativeRobustSumsFullRange->n - 1;i >= 0;i-- ) {
+            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." );
+        }
     /**************************************************************************
     Determine the mode in the range LQ:UQ.
-    **************************************************************************/
+    **************************************************************************/ 
     // Determine the bin with the peak value in the range LQ to UQ.
     maxBinNum = LQBinNum;
-    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
-    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];
-    }
-
+    maxBinCount = robustHistogramVector->data.F32[ LQBinNum ];
+    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 ];
+        }
+        
     // 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);
-
+    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
     **************************************************************************/
-    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    for (i=0;i<robustHistogramVector->n;i++) {
-        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
-    }
+    cumulativeRobustSumsDlRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 );
+    for ( i = 0;i < robustHistogramVector->n;i++ ) {
+            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];
-        }
-    }
-
+    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 ];
+                }
+        }
+        
     // 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
     // that bin (this is a non-exact approximation).
     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];
-        }
-    }
-    myMean/= countFloat;
-
+    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 ];
+                }
+        }
+    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];
-        }
-    }
-    myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1));
-
+    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 ];
+                }
+        }
+    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;
-    }
-
-    if  (stats->options & PS_STAT_ROBUST_MODE) {
-        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;
-    }
-
+    if ( stats->options & PS_STAT_ROBUST_MEAN ) {
+            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 ] );
+        }
+        
+    if ( stats->options & PS_STAT_ROBUST_STDEV ) {
+            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 +1374 @@
     // 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 psGetArrayPolynomial().  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 = psGetArrayPolynomial(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 ( 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 ] );
+                }
+        }
+        
     // The lower/upper quartile calculations are very similar to the median
     // calculations.  We fit a quadratic to the array containing the
@@ -1425 +1425 @@
     // f(x) equals the lower/upper quartile exactly.
     //
-    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 = psGetArrayPolynomial(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 = psGetArrayPolynomial(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]);
-        }
-    }
+    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 ] );
+                }
+        }
     stats->robustNfit = sumNfit;
     stats->robustN50 = sumN50;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(tmpStats);
-    psFree(robustHistogram);
-    psFree(myPoly);
-    psFree(cumulativeRobustSumsFullRange);
-    psFree(cumulativeRobustSumsDlRange);
+    
+    psFree( x );
+    psFree( y );
+    psFree( yErr );
+    psFree( tmpStats );
+    psFree( robustHistogram );
+    psFree( myPoly );
+    psFree( cumulativeRobustSumsFullRange );
+    psFree( cumulativeRobustSumsDlRange );
 }
 
@@ -1568 +1568 @@
 /*****************************************************************************/
 
-static void histogramFree(psHistogram *myHist);
+static void histogramFree( psHistogram *myHist );
 
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
  *****************************************************************************/
-psStats *psStatsAlloc(psStatsOptions options)
-{
-    psStats *newStruct = NULL;
-
-    newStruct = (psStats *) psAlloc(sizeof(psStats));
+psStats *psStatsAlloc( psStatsOptions options )
+{
+    psStats * newStruct = NULL;
+    
+    newStruct = ( psStats * ) psAlloc( sizeof( psStats ) );
     newStruct->sampleMean = NAN;
     newStruct->sampleMedian = NAN;
@@ -1600 +1600 @@
     newStruct->binsize = NAN;
     newStruct->options = options;
-
-    return(newStruct);
+    
+    return ( newStruct );
 }
 
@@ -1616 +1616 @@
     The histogram structure
  *****************************************************************************/
-psHistogram *psHistogramAlloc(float lower,
-                              float upper,
-                              int n)
+psHistogram *psHistogramAlloc( float lower,
+                               float upper,
+                               int n )
 {
     int i = 0;                          // Loop index variable
     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);
-    }
-    if (n < 0) {
-        psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n);
-    }
-
+    if ( n == 0 ) {
+            return ( NULL );
+        }
+    if ( n < 0 ) {
+            psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n );
+        }
+        
     // NOTE: Verify that this is the correct action.
-    if (lower > upper) {
-        return(NULL);
-    }
-
+    if ( lower > upper ) {
+            return ( NULL );
+        }
+        
     // Allocate memory for the new histogram structure.  If there are N
     // bins, then there are N+1 bounds to those bins.
-    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    p_psMemSetDeallocator(newHist,(psFreeFcn)histogramFree);
-    newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32);
+    newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) );
+    p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree );
+    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;
+    binSize = ( upper - lower ) / ( float ) n;
     // NOTE: Is the following necessary?  It prevents the max data point
     // from being in a non-existant bin.
-    binSize+= FLT_EPSILON;
-    for (i=0;i<n+1;i++) {
-        newHist->bounds->data.F32[i] = lower + (binSize * (float) i);
-    }
-
+    binSize += FLT_EPSILON;
+    for ( i = 0;i < n + 1;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 = 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;
-    }
-
+    for ( i = 0;i < newHist->nums->n;i++ ) {
+            newHist->nums->data.U32[ i ] = 0;
+        }
+        
     // Initialize the other members.
     newHist->minNum = 0;
     newHist->maxNum = 0;
     newHist->uniform = true;
-
-    return(newHist);
+    
+    return ( newHist );
 }
 
@@ -1677 +1677 @@
     The histogram structure
  *****************************************************************************/
-psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
-{
-    psHistogram *newHist = NULL;        // The new histogram structure
+psHistogram *psHistogramAllocGeneric( const psVector *restrict bounds )
+{
+    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);
-    }
-
+    if ( bounds == 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);
-    }
-
-    if (bounds->type.type != PS_TYPE_F32) {
-        // psAbort(__func__, "psHistogram request a bound which is not type F32");
-        return(NULL);
-    }
-
+    if ( bounds->n <= 1 ) {
+            // 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 );
+        }
+        
     // Allocate memory for the new histogram structure.
-    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    p_psMemSetDeallocator(newHist,(psFreeFcn)histogramFree);
-    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
+    newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) );
+    p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree );
+    newHist->bounds = psVectorAlloc( bounds->n, PS_TYPE_F32 );
     newHist->bounds->n = newHist->bounds->nalloc;
-    for (i=0;i<bounds->n;i++) {
-        newHist->bounds->data.F32[i] = bounds->data.F32[i];
-    }
-
+    for ( i = 0;i < bounds->n;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.
-    newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_U32);
+    newHist->nums = psVectorAlloc( ( bounds->n ) - 1, PS_TYPE_U32 );
     newHist->nums->n = newHist->nums->nalloc;
-    for (i=0;i<newHist->nums->n;i++) {
-        newHist->nums->data.U32[i] = 0;
-    }
-
+    for ( i = 0;i < newHist->nums->n;i++ ) {
+            newHist->nums->data.U32[ i ] = 0;
+        }
+        
     // Initialize the other members.
     newHist->minNum = 0;
     newHist->maxNum = 0;
     newHist->uniform = false;
-
-    return(newHist);
-}
-
-static void histogramFree(psHistogram *myHist)
-{
-    psFree(myHist->bounds);
-    psFree(myHist->nums);
+    
+    return ( newHist );
+}
+
+static void histogramFree( psHistogram *myHist )
+{
+    psFree( myHist->bounds );
+    psFree( myHist->nums );
 }
 
@@ -1745 +1745 @@
     The histogram structure "out".
  *****************************************************************************/
-psHistogram *psVectorHistogram(psHistogram *out,
-                               const psVector *restrict in,
-                               const psVector *restrict mask,
-                               unsigned int maskVal)
+psHistogram *psVectorHistogram( psHistogram *out,
+                                const psVector *restrict in,
+                                const psVector *restrict mask,
+                                unsigned int maskVal )
 {
     int i = 0;                                  // Loop index variable
@@ -1755 +1755 @@
     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);
-    }
-
+    if ( out == 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.");
-    }
-
-    if (out->nums->type.type != PS_TYPE_U32) {
-        psAbort(__func__,
-                "Only data type PS_TYPE_U32 for output.nums member.");
-    }
-
+    if ( out->bounds->type.type != PS_TYPE_F32 ) {
+            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." );
+        }
+        
     // NOTE: Verify that this is the correct action.
-    if (in == NULL) {
-        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__,
-                    "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");
-        }
-    }
+    if ( in == NULL ) {
+            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__,
+                             "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])++;
-                        }
-                    }
-                }
-            }
-        }
-    }
-    return(out);
+    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 ] ) ++;
+                                                }
+                                        }
+                                }
+                        }
+                }
+        }
+    return ( out );
 }
 
@@ -1852 +1852 @@
 the various stat functions.
  *****************************************************************************/
-psVector *p_psConvertToF32(psStats *stats,
-                           psVector *in,
-                           psVector *mask,
-                           unsigned int maskVal)
+psVector *p_psConvertToF32( psStats *stats,
+                            psVector *in,
+                            psVector *mask,
+                            unsigned int maskVal )
 {
     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.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);
-    }
-    return(tmp);
+    
+    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.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 );
+        }
+    return ( tmp );
 }
 
@@ -1901 +1901 @@
 macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
  *****************************************************************************/
-psStats *psVectorStats(psStats *stats,
-                       psVector *in,
-                       psVector *mask,
-                       unsigned int maskVal)
-{
-    psVector *inF32;
+psStats *psVectorStats( psStats *stats,
+                        psVector *in,
+                        psVector *mask,
+                        unsigned int maskVal )
+{
+    psVector * inF32;
     int mustFreeTmp = 1;
-
+    
     // NOTE: Verify that this is the correct action.
-    if (in == NULL) {
-        return(stats);
-    }
-    if (stats == NULL) {
-        return(NULL);
-    }
-
-    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
-    if (inF32 == NULL) {
-        inF32 = in;
-        mustFreeTmp = 0;
-    }
-
+    if ( in == NULL ) {
+            return ( stats );
+        }
+    if ( stats == NULL ) {
+            return ( NULL );
+        }
+        
+    inF32 = p_psConvertToF32( stats, in, mask, maskVal );
+    if ( inF32 == NULL ) {
+            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);
-    }
-
+    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 );
+        }
+        
     //    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);
-    }
-
+    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 );
+        }
+        
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_MEAN) {
-        p_psVectorSampleMean(in, mask, maskVal, stats);
-    }
-
+    if ( stats->options & PS_STAT_SAMPLE_MEAN ) {
+            p_psVectorSampleMean( in, mask, maskVal, stats );
+        }
+        
     // ************************************************************************
-    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
-        p_psVectorSampleMedian(in, mask, maskVal, stats);
-    }
-
+    if ( stats->options & PS_STAT_SAMPLE_MEDIAN ) {
+            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);
-    }
-
+    if ( stats->options & PS_STAT_SAMPLE_STDEV ) {
+            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);
-    }
-
+    if ( stats->options & PS_STAT_SAMPLE_QUARTILE ) {
+            p_psVectorSampleQuartiles( in, mask, maskVal, stats );
+        }
+        
     // Since the various robust stats quantities share much computation, they
     // are grouped together in a single private function:
     // p_psVectorRobustStats()
-    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
-            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
-            (stats->options & PS_STAT_ROBUST_MODE) ||
-            (stats->options & PS_STAT_ROBUST_STDEV) ||
-            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
-        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);
-    }
-
+    if ( ( stats->options & PS_STAT_ROBUST_MEAN ) ||
+            ( stats->options & PS_STAT_ROBUST_MEDIAN ) ||
+            ( stats->options & PS_STAT_ROBUST_MODE ) ||
+            ( stats->options & PS_STAT_ROBUST_STDEV ) ||
+            ( stats->options & PS_STAT_ROBUST_QUARTILE ) ) {
+            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 );
+        }
+        
     // ************************************************************************
-    if (stats->options & PS_STAT_MAX) {
-        p_psVectorMax(in, mask, maskVal, stats);
-    }
-
+    if ( stats->options & PS_STAT_MAX ) {
+            p_psVectorMax( in, mask, maskVal, stats );
+        }
+        
     // ************************************************************************
-    if (stats->options & PS_STAT_MIN) {
-        p_psVectorMin(in, mask, maskVal, stats);
-    }
-
-    if (mustFreeTmp == 1) {
-        psFree(inF32);
-    }
-    return(stats);
-}
+    if ( stats->options & PS_STAT_MIN ) {
+            p_psVectorMin( in, mask, maskVal, stats );
+        }
+        
+    if ( mustFreeTmp == 1 ) {
+            psFree( inF32 );
+        }
+    return ( stats );
+}