IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 6, 2004, 2:06:06 PM (22 years ago)
Author:
desonia
Message:

another attempt to get astyle to get it right.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psStats.c

    r1406 r1407  
     1
    12/** @file  psStats.c
    23 *  \brief basic statistical operations
    … …  
    910 *  @author George Gusciora, MHPCC
    1011 *
    11  *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-08-06 22:34:05 $
     12 *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-08-07 00:06:06 $
    1314 *
    1415 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2122#include <float.h>
    2223#include <math.h>
     24
    2325/*****************************************************************************/
     26
    2427/* INCLUDE FILES                                                             */
     28
    2529/*****************************************************************************/
    2630#include "psMemory.h"
    … …  
    3539
    3640/*****************************************************************************/
     41
    3742/* DEFINE STATEMENTS                                                         */
     43
    3844/*****************************************************************************/
    3945// will use robust statistical methods.
    40 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
     46#define GAUSS_WIDTH 20                     // The width of the Gaussian or boxcar smoothing.
    4147#define CLIPPED_NUM_ITER_LB 1
    4248#define CLIPPED_NUM_ITER_UB 10
    … …  
    4854#define MAX_ITERATIONS 10
    4955
    50 void p_psVectorRobustStats( const psVector *restrict myVector,
    51                             const psVector *restrict maskVector,
    52                             unsigned int maskVal,
    53                             psStats *stats );
     56void p_psVectorRobustStats(const psVector * restrict myVector,
     57                           const psVector * restrict maskVector, unsigned int maskVal, psStats * stats);
    5458
    5559/** Preprocessor macro to generate error on an incorrect type */
    … …  
    8387printf("\n"); \
    8488
     89
    8590/*****************************************************************************/
     91
    8692/* TYPE DEFINITIONS                                                          */
     93
    8794/*****************************************************************************/
    8895
    8996/*****************************************************************************/
     97
    9098/* GLOBAL VARIABLES                                                          */
     99
    91100/*****************************************************************************/
    92101
    … …  
    94103
    95104/*****************************************************************************/
     105
    96106/* FILE STATIC VARIABLES                                                     */
     107
    97108/*****************************************************************************/
    98109
    … …  
    100111
    101112/*****************************************************************************/
     113
    102114/* FUNCTION IMPLEMENTATION - LOCAL                                           */
     115
    103116/*****************************************************************************/
    104117
    105 bool p_psGetStatValue( const psStats* stats, double* value )
    106 {
    107 
    108     switch ( stats->options &
    109              ~ ( PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE ) ) {
     118bool p_psGetStatValue(const psStats * stats, double *value)
     119{
     120
     121    switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) {
    110122    case PS_STAT_SAMPLE_MEAN:
    111123        *value = stats->sampleMean;
    … …  
    163175 *****************************************************************************/
    164176
    165 void p_psVectorPrint( psVector *myVector,
    166                       psVector *maskVector,
    167                       unsigned int maskVal,
    168                       psStats *stats )
    169 {
    170     int i = 0;                                  // Loop index variable.
    171 
    172     for ( i = 0;i < myVector->n;i++ ) {
    173         if ( maskVector != NULL )
    174             printf( "Element %d is %f (mask is %d)\n", i, myVector->data.F32[ i ], maskVector->data.U8[ i ] );
     177void p_psVectorPrint(psVector * myVector, psVector * maskVector, unsigned int maskVal, psStats * stats)
     178{
     179    int i = 0;                  // Loop index variable.
     180
     181    for (i = 0; i < myVector->n; i++) {
     182        if (maskVector != NULL)
     183            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
    175184        else
    176             printf( "Element %d is %f\n", i, myVector->data.F32[ i ] );
     185            printf("Element %d is %f\n", i, myVector->data.F32[i]);
    177186    }
    178187}
    … …  
    195204 *****************************************************************************/
    196205
    197 
    198206/******************************************************************************
    199207p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
    … …  
    210218 *****************************************************************************/
    211219
    212 void p_psVectorSampleMean( const psVector *restrict myVector,
    213                            const psVector *restrict maskVector,
    214                            unsigned int maskVal,
    215                            psStats *stats )
    216 {
    217     int i = 0;                                  // Loop index variable
    218     float mean = 0.0;                           // The mean
    219     int count = 0;                              // # of points in this mean?
    220     float rangeMin = 0.0;                       // Exclude data below this
    221     float rangeMax = 0.0;                       // Exclude date above this
     220void p_psVectorSampleMean(const psVector * restrict myVector,
     221                          const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     222{
     223    int i = 0;                  // Loop index variable
     224    float mean = 0.0;           // The mean
     225    int count = 0;              // # of points in this mean?
     226    float rangeMin = 0.0;       // Exclude data below this
     227    float rangeMax = 0.0;       // Exclude date above this
    222228
    223229    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
    224230    // loop.
    225     if ( stats->options & PS_STAT_USE_RANGE ) {
     231    if (stats->options & PS_STAT_USE_RANGE) {
    226232        rangeMin = stats->min;
    227233        rangeMax = stats->max;
    228         if ( maskVector != NULL ) {
    229             for ( i = 0;i < myVector->n;i++ ) {
     234        if (maskVector != NULL) {
     235            for (i = 0; i < myVector->n; i++) {
    230236                // Check if the data is with the specified range
    231                 if ( !( maskVal & maskVector->data.U8[ i ] ) &&
    232                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    233                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    234                     mean += myVector->data.F32[ i ];
     237                if (!(maskVal & maskVector->data.U8[i]) &&
     238                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     239                    mean += myVector->data.F32[i];
    235240                    count++;
    236241                }
    237242            }
    238             mean /= ( float ) count;
     243            mean /= (float)count;
    239244        } else {
    240             for ( i = 0;i < myVector->n;i++ ) {
    241                 if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
    242                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    243                     mean += myVector->data.F32[ i ];
     245            for (i = 0; i < myVector->n; i++) {
     246                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     247                    mean += myVector->data.F32[i];
    244248                    count++;
    245249                }
    246250            }
    247             mean /= ( float ) count;
     251            mean /= (float)count;
    248252        }
    249253    } else {
    250         if ( maskVector != NULL ) {
    251             for ( i = 0;i < myVector->n;i++ ) {
    252                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    253                     mean += myVector->data.F32[ i ];
     254        if (maskVector != NULL) {
     255            for (i = 0; i < myVector->n; i++) {
     256                if (!(maskVal & maskVector->data.U8[i])) {
     257                    mean += myVector->data.F32[i];
    254258                    count++;
    255259                }
    256260            }
    257             mean /= ( float ) count;
     261            mean /= (float)count;
    258262        } else {
    259             for ( i = 0;i < myVector->n;i++ ) {
    260                 mean += myVector->data.F32[ i ];
    261             }
    262             mean /= ( float ) myVector->n;
     263            for (i = 0; i < myVector->n; i++) {
     264                mean += myVector->data.F32[i];
     265            }
     266            mean /= (float)myVector->n;
    263267        }
    264268    }
    … …  
    278282    NULL
    279283 *****************************************************************************/
    280 void p_psVectorMax( const psVector *restrict myVector,
    281                     const psVector *restrict maskVector,
    282                     unsigned int maskVal,
    283                     psStats *stats )
    284 {
    285     int i = 0;                                  // Loop index variable
    286     float max = -MY_MAX_FLOAT;                    // The calculated maximum
    287     float rangeMin = 0.0;                       // Exclude data below this
    288     float rangeMax = 0.0;                       // Exclude date above this
     284void p_psVectorMax(const psVector * restrict myVector,
     285                   const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     286{
     287    int i = 0;                  // Loop index variable
     288    float max = -MY_MAX_FLOAT;  // The calculated maximum
     289    float rangeMin = 0.0;       // Exclude data below this
     290    float rangeMax = 0.0;       // Exclude date above this
    289291
    290292    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    291     if ( stats->options & PS_STAT_USE_RANGE ) {
     293    if (stats->options & PS_STAT_USE_RANGE) {
    292294        rangeMin = stats->min;
    293295        rangeMax = stats->max;
    294         if ( maskVector != NULL ) {
    295             for ( i = 0;i < myVector->n;i++ ) {
    296                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    297                     if ( ( myVector->data.F32[ i ] > max ) &&
    298                             ( rangeMin <= myVector->data.F32[ i ] ) &&
    299                             ( myVector->data.F32[ i ] <= rangeMax ) ) {
    300                         max = myVector->data.F32[ i ];
     296        if (maskVector != NULL) {
     297            for (i = 0; i < myVector->n; i++) {
     298                if (!(maskVal & maskVector->data.U8[i])) {
     299                    if ((myVector->data.F32[i] > max) &&
     300                            (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     301                        max = myVector->data.F32[i];
    301302                    }
    302303                }
    303304            }
    304305        } else {
    305             for ( i = 0;i < myVector->n;i++ ) {
    306                 if ( ( myVector->data.F32[ i ] > max ) &&
    307                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    308                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    309                     max = myVector->data.F32[ i ];
     306            for (i = 0; i < myVector->n; i++) {
     307                if ((myVector->data.F32[i] > max) &&
     308                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     309                    max = myVector->data.F32[i];
    310310                }
    311311            }
    312312        }
    313313    } else {
    314         if ( maskVector != NULL ) {
    315             for ( i = 0;i < myVector->n;i++ ) {
    316                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    317                     if ( myVector->data.F32[ i ] > max ) {
    318                         max = myVector->data.F32[ i ];
     314        if (maskVector != NULL) {
     315            for (i = 0; i < myVector->n; i++) {
     316                if (!(maskVal & maskVector->data.U8[i])) {
     317                    if (myVector->data.F32[i] > max) {
     318                        max = myVector->data.F32[i];
    319319                    }
    320320                }
    321321            }
    322322        } else {
    323             for ( i = 0;i < myVector->n;i++ ) {
    324                 if ( myVector->data.F32[ i ] > max ) {
    325                     max = myVector->data.F32[ i ];
     323            for (i = 0; i < myVector->n; i++) {
     324                if (myVector->data.F32[i] > max) {
     325                    max = myVector->data.F32[i];
    326326                }
    327327            }
    … …  
    343343    NULL
    344344 *****************************************************************************/
    345 void p_psVectorMin( const psVector *restrict myVector,
    346                     const psVector *restrict maskVector,
    347                     unsigned int maskVal,
    348                     psStats *stats )
    349 {
    350     int i = 0;                                  // Loop index variable
    351     float min = MY_MAX_FLOAT;                     // The calculated maximum
    352     float rangeMin = 0.0;                       // Exclude data below this
    353     float rangeMax = 0.0;                       // Exclude date above this
    354 
    355     if ( stats->options & PS_STAT_USE_RANGE ) {
     345void p_psVectorMin(const psVector * restrict myVector,
     346                   const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     347{
     348    int i = 0;                  // Loop index variable
     349    float min = MY_MAX_FLOAT;   // The calculated maximum
     350    float rangeMin = 0.0;       // Exclude data below this
     351    float rangeMax = 0.0;       // Exclude date above this
     352
     353    if (stats->options & PS_STAT_USE_RANGE) {
    356354        rangeMin = stats->min;
    357355        rangeMax = stats->max;
    358         if ( maskVector != NULL ) {
    359             for ( i = 0;i < myVector->n;i++ ) {
    360                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    361                     if ( ( myVector->data.F32[ i ] < min ) &&
    362                             ( rangeMin <= myVector->data.F32[ i ] ) &&
    363                             ( myVector->data.F32[ i ] <= rangeMax ) ) {
    364                         min = myVector->data.F32[ i ];
     356        if (maskVector != NULL) {
     357            for (i = 0; i < myVector->n; i++) {
     358                if (!(maskVal & maskVector->data.U8[i])) {
     359                    if ((myVector->data.F32[i] < min) &&
     360                            (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     361                        min = myVector->data.F32[i];
    365362                    }
    366363                }
    367364            }
    368365        } else {
    369             for ( i = 0;i < myVector->n;i++ ) {
    370                 if ( ( myVector->data.F32[ i ] < min ) &&
    371                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    372                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    373                     min = myVector->data.F32[ i ];
     366            for (i = 0; i < myVector->n; i++) {
     367                if ((myVector->data.F32[i] < min) &&
     368                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     369                    min = myVector->data.F32[i];
    374370                }
    375371            }
    376372        }
    377373    } else {
    378         if ( maskVector != NULL ) {
    379             for ( i = 0;i < myVector->n;i++ ) {
    380                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    381                     if ( myVector->data.F32[ i ] < min ) {
    382                         min = myVector->data.F32[ i ];
     374        if (maskVector != NULL) {
     375            for (i = 0; i < myVector->n; i++) {
     376                if (!(maskVal & maskVector->data.U8[i])) {
     377                    if (myVector->data.F32[i] < min) {
     378                        min = myVector->data.F32[i];
    383379                    }
    384380                }
    385381            }
    386382        } else {
    387             for ( i = 0;i < myVector->n;i++ ) {
    388                 if ( myVector->data.F32[ i ] < min ) {
    389                     min = myVector->data.F32[ i ];
     383            for (i = 0; i < myVector->n; i++) {
     384                if (myVector->data.F32[i] < min) {
     385                    min = myVector->data.F32[i];
    390386                }
    391387            }
    … …  
    408404    NULL
    409405 *****************************************************************************/
    410 int p_psVectorNValues( const psVector *restrict myVector,
    411                        const psVector *restrict maskVector,
    412                        unsigned int maskVal,
    413                        psStats *stats )
    414 {
    415     int i = 0;                                  // Loop index variable
    416     int numData = 0;                            // The number of data points
    417     float rangeMin = 0.0;                       // Exclude data below this
    418     float rangeMax = 0.0;                       // Exclude date above this
    419 
    420     if ( stats->options & PS_STAT_USE_RANGE ) {
     406int p_psVectorNValues(const psVector * restrict myVector,
     407                      const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     408{
     409    int i = 0;                  // Loop index variable
     410    int numData = 0;            // The number of data points
     411    float rangeMin = 0.0;       // Exclude data below this
     412    float rangeMax = 0.0;       // Exclude date above this
     413
     414    if (stats->options & PS_STAT_USE_RANGE) {
    421415        rangeMin = stats->min;
    422416        rangeMax = stats->max;
    423         if ( maskVector != NULL ) {
    424             for ( i = 0;i < myVector->n;i++ ) {
    425                 if ( !( maskVal & maskVector->data.U8[ i ] ) &&
    426                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    427                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
     417        if (maskVector != NULL) {
     418            for (i = 0; i < myVector->n; i++) {
     419                if (!(maskVal & maskVector->data.U8[i]) &&
     420                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
    428421                    numData++;
    429422                }
    430423            }
    431424        } else {
    432             for ( i = 0;i < myVector->n;i++ ) {
    433                 if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
    434                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
     425            for (i = 0; i < myVector->n; i++) {
     426                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
    435427                    numData++;
    436428                }
    … …  
    440432        rangeMin = stats->min;
    441433        rangeMax = stats->max;
    442         if ( maskVector != NULL ) {
    443             for ( i = 0;i < myVector->n;i++ ) {
    444                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
     434        if (maskVector != NULL) {
     435            for (i = 0; i < myVector->n; i++) {
     436                if (!(maskVal & maskVector->data.U8[i])) {
    445437                    numData++;
    446438                }
    … …  
    450442        }
    451443    }
    452     return ( numData );
    453 }
    454 
    455 
     444    return (numData);
     445}
    456446
    457447/******************************************************************************
    … …  
    466456    NULL
    467457 *****************************************************************************/
    468 void p_psVectorSampleMedian( const psVector *restrict myVector,
    469                              const psVector *restrict maskVector,
    470                              unsigned int maskVal,
    471                              psStats *stats )
    472 {
    473     psVector * unsortedVector = NULL;            // Temporary vector
    474     psVector *sortedVector = NULL;              // Temporary vector
    475     int i = 0;                                  // Loop index variable
    476     int count = 0;                              // # of points in this mean?
    477     int nValues = 0;                            // # of points in vector
    478     float rangeMin = 0.0;                       // Exclude data below this
    479     float rangeMax = 0.0;                       // Exclude date above this
    480 
     458void p_psVectorSampleMedian(const psVector * restrict myVector,
     459                            const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     460{
     461    psVector *unsortedVector = NULL;    // Temporary vector
     462    psVector *sortedVector = NULL;      // Temporary vector
     463    int i = 0;                  // Loop index variable
     464    int count = 0;              // # of points in this mean?
     465    int nValues = 0;            // # of points in vector
     466    float rangeMin = 0.0;       // Exclude data below this
     467    float rangeMax = 0.0;       // Exclude date above this
    481468
    482469    // Determine if the number of data points exceed a threshold which will
    … …  
    485472    // regardless of the vector size.
    486473    /*
    487         if (myVector->n > stats->sampleLimit) {
    488             psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    489      
    490             // Calculate the robust quartiles.
    491             stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    492             p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    493      
    494             // Store the robust quartiles into the sample quartile members.
    495             stats->sampleMedian = stats2->robustMedian;
    496      
    497             // Free temporary data buffers.
    498             psFree(stats2);
    499      
    500             // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    501             stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    502      
    503             return;
    504         }
    505     */
     474     * if (myVector->n > stats->sampleLimit) { psAbort(__func__, "Robust Statistic Algorithms have not yet
     475     * been defined or implemented.");
     476     *
     477     * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     478     * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     479     *
     480     * // Store the robust quartiles into the sample quartile members. stats->sampleMedian =
     481     * stats2->robustMedian;
     482     *
     483     * // Free temporary data buffers. psFree(stats2);
     484     *
     485     * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options |
     486     * PS_STAT_ROBUST_FOR_SAMPLE;
     487     *
     488     * return; } */
    506489
    507490    // Determine how many data points fit inside this min/max range
    508491    // and are not masked, IF the maskVector is not NULL>
    509     nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats );
     492    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    510493
    511494    // Allocate temporary vectors for the data.
    512     unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     495    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    513496    unsortedVector->n = unsortedVector->nalloc;
    514497
    515     sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     498    sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    516499    sortedVector->n = sortedVector->nalloc;
    517500
    518501    // Determine if we must only use data points within a min/max range.
    519     if ( stats->options & PS_STAT_USE_RANGE ) {
     502    if (stats->options & PS_STAT_USE_RANGE) {
    520503        rangeMin = stats->min;
    521504        rangeMax = stats->max;
    … …  
    524507        // into the temporary vectors.
    525508        count = 0;
    526         if ( maskVector != NULL ) {
    527             for ( i = 0;i < myVector->n;i++ ) {
    528                 if ( !( maskVal & maskVector->data.U8[ i ] ) &&
    529                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    530                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    531                     unsortedVector->data.F32[ count++ ] = maskVector->data.F32[ i ];
     509        if (maskVector != NULL) {
     510            for (i = 0; i < myVector->n; i++) {
     511                if (!(maskVal & maskVector->data.U8[i]) &&
     512                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     513                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
    532514                }
    533515            }
    534516        } else {
    535             for ( i = 0;i < myVector->n;i++ ) {
    536                 if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
    537                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    538                     unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
     517            for (i = 0; i < myVector->n; i++) {
     518                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     519                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
    539520                }
    540521            }
    … …  
    543524        // Store all non-masked data points into the temporary vectors.
    544525        count = 0;
    545         if ( maskVector != NULL ) {
    546             for ( i = 0;i < myVector->n;i++ ) {
    547                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    548                     unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
     526        if (maskVector != NULL) {
     527            for (i = 0; i < myVector->n; i++) {
     528                if (!(maskVal & maskVector->data.U8[i])) {
     529                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
    549530                }
    550531            }
    551532        } else {
    552             for ( i = 0;i < myVector->n;i++ ) {
    553                 unsortedVector->data.F32[ i ] = myVector->data.F32[ i ];
     533            for (i = 0; i < myVector->n; i++) {
     534                unsortedVector->data.F32[i] = myVector->data.F32[i];
    554535            }
    555536        }
    556537    }
    557538    // Sort the temporary vectors.
    558     psVectorSort( sortedVector, unsortedVector );
     539    psVectorSort(sortedVector, unsortedVector);
    559540
    560541    // Calculate the median exactly.
    561542    // XXX: Is this the correct action?
    562     if ( 0 == ( nValues % 2 ) ) {
    563         stats->sampleMedian = 0.5 * ( sortedVector->data.F32[ ( nValues / 2 ) - 1 ] +
    564                                       sortedVector->data.F32[ nValues / 2 ] );
     543    if (0 == (nValues % 2)) {
     544        stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] +
     545                                     sortedVector->data.F32[nValues / 2]);
    565546    } else {
    566         stats->sampleMedian = sortedVector->data.F32[ nValues / 2 ];
     547        stats->sampleMedian = sortedVector->data.F32[nValues / 2];
    567548    }
    568549
    569550    // Free the temporary data structures.
    570     psFree( unsortedVector );
    571     psFree( sortedVector );
     551    psFree(unsortedVector);
     552    psFree(sortedVector);
    572553}
    573554
    … …  
    584565    XXX: use a static variable for gaussianCoefs[] and compute them once.
    585566 *****************************************************************************/
    586 psVector *p_psVectorsmoothHistGaussian( psHistogram *robustHistogram,
    587                                         float sigma )
    588 {
    589     int i = 0;                                  // Loop index variable
    590     int j = 0;                                  // Loop index variable
    591     float denom = 0.0;                           // Temporary variable
    592     float expo = 0.0;                            // Temporary variable
    593     float gaussianCoefs[ 1 + ( 2 * GAUSS_WIDTH ) ]; // The Gaussian Coefficients
    594     psVector *smooth = psVectorAlloc( robustHistogram->nums->n, PS_TYPE_F32 );
    595 
    596     for ( i = 0;i < ( 1 + ( 2 * GAUSS_WIDTH ) );i++ ) {
    597         if ( fabs( sigma ) >= FLT_EPSILON ) {
     567psVector *p_psVectorsmoothHistGaussian(psHistogram * robustHistogram, float sigma)
     568{
     569    int i = 0;                  // Loop index variable
     570    int j = 0;                  // Loop index variable
     571    float denom = 0.0;          // Temporary variable
     572    float expo = 0.0;           // Temporary variable
     573    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
     574    psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
     575
     576    for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) {
     577        if (fabs(sigma) >= FLT_EPSILON) {
    598578            // If sigma does not equal zero, then we use Gaussian smoothing.
    599579            #ifdef  DARWIN
    600             denom = ( float ) sqrt( 2.0 * M_PI * sigma * sigma );
     580            denom = (float)sqrt(2.0 * M_PI * sigma * sigma);
    601581            #else
    602582
    603             denom = sqrtf( 2.0 * M_PI * sigma * sigma );
     583            denom = sqrtf(2.0 * M_PI * sigma * sigma);
    604584            #endif
    605585
    606             expo = - ( float ) ( ( i - GAUSS_WIDTH ) * ( i - GAUSS_WIDTH ) );
    607             expo /= ( 2.0 * sigma * sigma );
    608             gaussianCoefs[ i ] = exp( expo / denom );
     586            expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH));
     587            expo /= (2.0 * sigma * sigma);
     588            gaussianCoefs[i] = exp(expo / denom);
    609589
    610590            // NOTE: Gaussian smoothing just isn't working with low sigma
    … …  
    612592            // all zero, except for the middle coefficient, which is exactly
    613593            // one.  Therefore, I'm using boxcar smoothing.
    614             gaussianCoefs[ i ] = 1.0 / ( 1.0 + ( 2.0 * ( float ) GAUSS_WIDTH ) );
    615             //            printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
     594            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH));
     595            // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
    616596        } else {
    617             /* If sigma equals zero (all pixels have the same value)
    618              * the above code will divide by zero.  Therefore, we don't need
    619              * to smooth the data.
    620              */
    621             for ( i = 0;i < robustHistogram->nums->n;i++ ) {
    622                 smooth->data.F32[ i ] = ( float ) robustHistogram->nums->data.S32[ i ];
    623             }
    624             return ( smooth );
     597            /* If sigma equals zero (all pixels have the same value) the above code will divide by zero.
     598             * Therefore, we don't need to smooth the data. */
     599            for (i = 0; i < robustHistogram->nums->n; i++) {
     600                smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i];
     601            }
     602            return (smooth);
    625603        }
    626604    }
    627605
    628606    // Perform the actual smoothing.
    629     for ( i = 0;i < robustHistogram->nums->n;i++ ) {
    630         smooth->data.F32[ i ] = 0.0;
    631         for ( j = -GAUSS_WIDTH;j <= + GAUSS_WIDTH;j++ ) {
    632             if ( ( ( j + i ) >= 0 ) && ( ( j + i ) < smooth->n ) ) {
    633                 smooth->data.F32[ i ] += ( gaussianCoefs[ j + GAUSS_WIDTH ] *
    634                                            ( float ) robustHistogram->nums->data.S32[ j + i ] );
    635             }
    636         }
    637     }
    638     return ( smooth );
     607    for (i = 0; i < robustHistogram->nums->n; i++) {
     608        smooth->data.F32[i] = 0.0;
     609        for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) {
     610            if (((j + i) >= 0) && ((j + i) < smooth->n)) {
     611                smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] *
     612                                        (float)robustHistogram->nums->data.S32[j + i]);
     613            }
     614        }
     615    }
     616    return (smooth);
    639617}
    640618
    … …  
    650628    NULL
    651629 *****************************************************************************/
    652 void p_psVectorSampleQuartiles( const psVector *restrict myVector,
    653                                 const psVector *restrict maskVector,
    654                                 unsigned int maskVal,
    655                                 psStats *stats )
    656 {
    657     psVector * unsortedVector = NULL;            // Temporary vector
    658     psVector *sortedVector = NULL;              // Temporary vector
    659     int i = 0;                                  // Loop index variable
    660     int count = 0;                              // # of points in this mean?
    661     int nValues = 0;                            // # data points
    662     float rangeMin = 0.0;                       // Exclude data below this
    663     float rangeMax = 0.0;                       // Exclude date above this
     630void p_psVectorSampleQuartiles(const psVector * restrict myVector,
     631                               const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     632{
     633    psVector *unsortedVector = NULL;    // Temporary vector
     634    psVector *sortedVector = NULL;      // Temporary vector
     635    int i = 0;                  // Loop index variable
     636    int count = 0;              // # of points in this mean?
     637    int nValues = 0;            // # data points
     638    float rangeMin = 0.0;       // Exclude data below this
     639    float rangeMax = 0.0;       // Exclude date above this
    664640
    665641    // Determine how many data points fit inside this min/max range
    666642    // and are not maxed, IF the maskVector is not NULL>
    667     nValues = p_psVectorNValues( myVector, maskVector, maskVal, stats );
     643    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    668644
    669645    // Allocate temporary vectors for the data.
    670     unsortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     646    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    671647    unsortedVector->n = unsortedVector->nalloc;
    672     sortedVector = psVectorAlloc( nValues, PS_TYPE_F32 );
     648    sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    673649    sortedVector->n = sortedVector->nalloc;
    674650
    675651    // Determine if we must only use data points within a min/max range.
    676     if ( stats->options & PS_STAT_USE_RANGE ) {
     652    if (stats->options & PS_STAT_USE_RANGE) {
    677653        rangeMin = stats->min;
    678654        rangeMax = stats->max;
    … …  
    680656        // into the temporary vectors.
    681657        count = 0;
    682         if ( maskVector != NULL ) {
    683             for ( i = 0;i < myVector->n;i++ ) {
    684                 if ( !( maskVal & maskVector->data.U8[ i ] ) &&
    685                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    686                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    687                     unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
     658        if (maskVector != NULL) {
     659            for (i = 0; i < myVector->n; i++) {
     660                if (!(maskVal & maskVector->data.U8[i]) &&
     661                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     662                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
    688663                }
    689664            }
    690665        } else {
    691             for ( i = 0;i < myVector->n;i++ ) {
    692                 if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
    693                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    694                     unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
     666            for (i = 0; i < myVector->n; i++) {
     667                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     668                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
    695669                }
    696670            }
    … …  
    699673        // Store all non-masked data points into the temporary vectors.
    700674        count = 0;
    701         if ( maskVector != NULL ) {
    702             for ( i = 0;i < myVector->n;i++ ) {
    703                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    704                     unsortedVector->data.F32[ count++ ] = myVector->data.F32[ i ];
     675        if (maskVector != NULL) {
     676            for (i = 0; i < myVector->n; i++) {
     677                if (!(maskVal & maskVector->data.U8[i])) {
     678                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
    705679                }
    706680            }
    707681        } else {
    708             for ( i = 0;i < myVector->n;i++ ) {
    709                 unsortedVector->data.F32[ i ] = myVector->data.F32[ i ];
     682            for (i = 0; i < myVector->n; i++) {
     683                unsortedVector->data.F32[i] = myVector->data.F32[i];
    710684            }
    711685        }
    … …  
    713687
    714688    // Sort the temporary vectors.
    715     psVectorSort( sortedVector, unsortedVector );
     689    psVectorSort(sortedVector, unsortedVector);
    716690
    717691    // Calculate the quartile points exactly.
    718     stats->sampleUQ = sortedVector->data.F32[ 3 * ( nValues / 4 ) ];
    719     stats->sampleLQ = sortedVector->data.F32[ nValues / 4 ];
     692    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     693    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    720694
    721695    // Free the temporary data structures.
    722     psFree( unsortedVector );
    723     psFree( sortedVector );
    724 }
    725 
     696    psFree(unsortedVector);
     697    psFree(sortedVector);
     698}
    726699
    727700/******************************************************************************
    … …  
    737710 
    738711 *****************************************************************************/
    739 void p_psVectorSampleStdev( const psVector *restrict myVector,
    740                             const psVector *restrict maskVector,
    741                             unsigned int maskVal,
    742                             psStats *stats )
    743 {
    744     int i = 0;                                  // Loop index variable
    745     int countInt = 0;                           // # of data points being used
    746     float countFloat = 0.0;                     // # of data points being used
    747     float mean = 0.0;                           // The mean
    748     float diff = 0.0;                           // Used in calculating stdev
    749     float sumSquares = 0.0;                     // temporary variable
    750     float sumDiffs = 0.0;                       // temporary variable
    751     float rangeMin = 0.0;                       // Exclude data below this
    752     float rangeMax = 0.0;                       // Exclude date above this
     712void p_psVectorSampleStdev(const psVector * restrict myVector,
     713                           const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     714{
     715    int i = 0;                  // Loop index variable
     716    int countInt = 0;           // # of data points being used
     717    float countFloat = 0.0;     // # of data points being used
     718    float mean = 0.0;           // The mean
     719    float diff = 0.0;           // Used in calculating stdev
     720    float sumSquares = 0.0;     // temporary variable
     721    float sumDiffs = 0.0;       // temporary variable
     722    float rangeMin = 0.0;       // Exclude data below this
     723    float rangeMax = 0.0;       // Exclude date above this
    753724
    754725    // This procedure requires the mean.  If it has not been already
    755726    // calculated, then call p_psVectorSampleMean()
    756     if ( 0 != isnan( stats->sampleMean ) ) {
    757         p_psVectorSampleMean( myVector, maskVector, maskVal, stats );
     727    if (0 != isnan(stats->sampleMean)) {
     728        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    758729    }
    759730    mean = stats->sampleMean;
    760731
    761     if ( stats->options & PS_STAT_USE_RANGE ) {
    762         if ( maskVector != NULL ) {
    763             for ( i = 0;i < myVector->n;i++ ) {
    764                 if ( !( maskVal & maskVector->data.U8[ i ] ) &&
    765                         ( rangeMin <= myVector->data.F32[ i ] ) &&
    766                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    767                     diff = myVector->data.F32[ i ] - mean;
    768                     sumSquares += ( diff * diff );
     732    if (stats->options & PS_STAT_USE_RANGE) {
     733        if (maskVector != NULL) {
     734            for (i = 0; i < myVector->n; i++) {
     735                if (!(maskVal & maskVector->data.U8[i]) &&
     736                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     737                    diff = myVector->data.F32[i] - mean;
     738                    sumSquares += (diff * diff);
    769739                    sumDiffs += diff;
    770740                    countInt++;
    … …  
    772742            }
    773743        } else {
    774             for ( i = 0;i < myVector->n;i++ ) {
    775                 if ( ( rangeMin <= myVector->data.F32[ i ] ) &&
    776                         ( myVector->data.F32[ i ] <= rangeMax ) ) {
    777                     diff = myVector->data.F32[ i ] - mean;
    778                     sumSquares += ( diff * diff );
     744            for (i = 0; i < myVector->n; i++) {
     745                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
     746                    diff = myVector->data.F32[i] - mean;
     747                    sumSquares += (diff * diff);
    779748                    sumDiffs += diff;
    780749                    countInt++;
    … …  
    784753        }
    785754    } else {
    786         if ( maskVector != NULL ) {
    787             for ( i = 0;i < myVector->n;i++ ) {
    788                 if ( !( maskVal & maskVector->data.U8[ i ] ) ) {
    789                     diff = myVector->data.F32[ i ] - mean;
    790                     sumSquares += ( diff * diff );
     755        if (maskVector != NULL) {
     756            for (i = 0; i < myVector->n; i++) {
     757                if (!(maskVal & maskVector->data.U8[i])) {
     758                    diff = myVector->data.F32[i] - mean;
     759                    sumSquares += (diff * diff);
    791760                    sumDiffs += diff;
    792761                    countInt++;
    … …  
    794763            }
    795764        } else {
    796             for ( i = 0;i < myVector->n;i++ ) {
    797                 diff = myVector->data.F32[ i ] - mean;
    798                 sumSquares += ( diff * diff );
     765            for (i = 0; i < myVector->n; i++) {
     766                diff = myVector->data.F32[i] - mean;
     767                sumSquares += (diff * diff);
    799768                sumDiffs += diff;
    800769                countInt++;
    … …  
    803772        }
    804773    }
    805     countFloat = ( float ) countInt;
     774    countFloat = (float)countInt;
    806775
    807776    #ifdef DARWIN
    808777
    809     stats->sampleStdev = ( float ) sqrt( ( sumSquares - ( sumDiffs *
    810                                            sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     778    stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    811779    #else
    812780
    813     stats->sampleStdev = sqrtf( ( sumSquares - ( sumDiffs *
    814                                   sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     781    stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    815782    #endif
    816783}
    … …  
    828795    NULL
    829796 *****************************************************************************/
    830 void p_psVectorClippedStats( const psVector *restrict myVector,
    831                              const psVector *restrict maskVector,
    832                              unsigned int maskVal,
    833                              psStats *stats )
    834 {
    835     int i = 0;                                  // Loop index variable
    836     int j = 0;                                  // Loop index variable
    837     float clippedMean = 0.0;                    // self-explanatory
    838     float clippedStdev = 0.0;                   // self-explanatory
    839     float oldStanMean = 0.0;                    // Temporary variable
    840     float oldStanStdev = 0.0;                   // Temporary variable
    841     psVector *tmpMask = NULL;                   // Temporary vector
     797void p_psVectorClippedStats(const psVector * restrict myVector,
     798                            const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     799{
     800    int i = 0;                  // Loop index variable
     801    int j = 0;                  // Loop index variable
     802    float clippedMean = 0.0;    // self-explanatory
     803    float clippedStdev = 0.0;   // self-explanatory
     804    float oldStanMean = 0.0;    // Temporary variable
     805    float oldStanStdev = 0.0;   // Temporary variable
     806    psVector *tmpMask = NULL;   // Temporary vector
    842807
    843808    // Endure that stats->clipIter is within the proper range.
    844     if ( !( ( CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
    845             ( stats->clipIter <= CLIPPED_NUM_ITER_UB ) ) ) {
    846         psAbort( __func__, "Unallowed value for clipIter (%d).\n",
    847                  stats->clipIter );
    848     }
    849 
     809    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
     810        psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
     811    }
    850812    // Endure that stats->clipSigma is within the proper range.
    851     if ( !( ( CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
    852             ( stats->clipSigma <= CLIPPED_SIGMA_UB ) ) ) {
    853         psAbort( __func__, "Unallowed value for clipSigma (%f).\n",
    854                  stats->clipSigma );
    855     }
    856 
     813    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
     814        psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
     815    }
    857816    // We allocate a temporary mask vector since during the iterative
    858817    // steps that follow, we will be masking off additional data points.
    859818    // However, we do no want to modify the original mask vector.
    860     tmpMask = psVectorAlloc( myVector->n, PS_TYPE_U8 );
     819    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    861820    tmpMask->n = myVector->n;
    862821
    863822    // If we were called with a mask vector, then initialize the temporary
    864823    // mask vector with those values.
    865     if ( maskVector != NULL ) {
    866         for ( i = 0;i < tmpMask->n;i++ ) {
    867             tmpMask->data.U8[ i ] = maskVector->data.U8[ i ];
    868         }
    869     }
    870 
     824    if (maskVector != NULL) {
     825        for (i = 0; i < tmpMask->n; i++) {
     826            tmpMask->data.U8[i] = maskVector->data.U8[i];
     827        }
     828    }
    871829    // 1. Compute the sample median.
    872830    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
    873831    // median here, not the mean.
    874     p_psVectorSampleMedian( myVector, maskVector, maskVal, stats );
     832    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
    875833
    876834    // 2. Compute the sample standard deviation.
    877     p_psVectorSampleStdev( myVector, maskVector, maskVal, stats );
     835    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
    878836
    879837    // 3. Use the sample median as the first estimator of the mean X.
    … …  
    889847
    890848    // 5. Repeat N times:
    891     for ( i = 0;i < stats->clipIter;i++ ) {
    892         for ( j = 0;j < myVector->n;j++ ) {
     849    for (i = 0; i < stats->clipIter; i++) {
     850        for (j = 0; j < myVector->n; j++) {
    893851            // a) Exclude all values x_i for which |x_i - x| > K * stdev
    894             if ( fabs( myVector->data.F32[ j ] - clippedMean ) >
    895                     ( stats->clipSigma * clippedStdev ) ) {
    896                 tmpMask->data.U8[ i ] = 0xff;
     852            if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
     853                tmpMask->data.U8[i] = 0xff;
    897854            }
    898855            // b) compute new mean and stdev
    899             p_psVectorSampleMedian( myVector, tmpMask, maskVal, stats );
    900             p_psVectorSampleStdev( myVector, tmpMask, maskVal, stats );
     856            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
     857            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
    901858
    902859            // c) Use the new mean for x
    … …  
    912869
    913870    // 7. The last calcuated value of x is the cliped mean.
    914     if ( stats->options & PS_STAT_CLIPPED_MEAN ) {
     871    if (stats->options & PS_STAT_CLIPPED_MEAN) {
    915872        stats->clippedMean = clippedMean;
    916873    }
    917 
    918874    // 8. The last calcuated value of stdev is the cliped stdev.
    919     if ( stats->options & PS_STAT_CLIPPED_STDEV ) {
     875    if (stats->options & PS_STAT_CLIPPED_STDEV) {
    920876        stats->clippedStdev = clippedStdev;
    921877    }
    922878
    923     psFree( tmpMask );
     879    psFree(tmpMask);
    924880}
    925881
    … …  
    928884elements of a vector to a range between 0.0 and 1.0.
    929885 *****************************************************************************/
    930 void p_psNormalizeVector( psVector *myData )
    931 {
    932     float min = ( float ) HUGE;
    933     float max = ( float ) - HUGE;
     886void p_psNormalizeVector(psVector * myData)
     887{
     888    float min = (float)HUGE;
     889    float max = (float)-HUGE;
    934890    float range = 0.0;
    935891    int i = 0;
    936892
    937     for ( i = 0;i < myData->n;i++ ) {
    938         if ( myData->data.F32[ i ] < min ) {
    939             min = myData->data.F32[ i ];
    940         }
    941         if ( myData->data.F32[ i ] > max ) {
    942             max = myData->data.F32[ i ];
     893    for (i = 0; i < myData->n; i++) {
     894        if (myData->data.F32[i] < min) {
     895            min = myData->data.F32[i];
     896        }
     897        if (myData->data.F32[i] > max) {
     898            max = myData->data.F32[i];
    943899        }
    944900    }
    945901
    946902    range = max - min;
    947     for ( i = 0;i < myData->n;i++ ) {
    948         myData->data.F32[ i ] = ( myData->data.F32[ i ] - min ) / range;
    949     }
    950 }
    951 
     903    for (i = 0; i < myData->n; i++) {
     904        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
     905    }
     906}
    952907
    953908/*****************************************************************************
    … …  
    956911specified data point.
    957912 *****************************************************************************/
    958 float p_psGaussian( const psVector *restrict myData,
    959                     const psVector *restrict myParams )
    960 {
    961     float x = myData->data.F32[ 0 ];
    962     float mean = myParams->data.F32[ 0 ];
    963     float stdev = myParams->data.F32[ 1 ];
    964     float tmp = exp( -( ( x - mean ) * ( x - mean ) ) / ( 2.0 * stdev * stdev ) );
    965     tmp /= ( ( float ) sqrt( 2.0 * M_PI * ( stdev * stdev ) ) );
    966 
    967     //    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
    968     return ( tmp );
     913float p_psGaussian(const psVector * restrict myData, const psVector * restrict myParams)
     914{
     915    float x = myData->data.F32[0];
     916    float mean = myParams->data.F32[0];
     917    float stdev = myParams->data.F32[1];
     918    float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev));
     919
     920    tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
     921
     922    // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     923    return (tmp);
    969924}
    970925
    … …  
    973928calculates the specified partial derivative of the above Gaussian function.
    974929 *****************************************************************************/
    975 float p_psGaussianDeriv( const psVector *restrict myData,
    976                          const psVector *restrict myParams,
    977                          int whichParam )
    978 {
    979     float x = myData->data.F32[ 0 ];
    980     float mean = myParams->data.F32[ 0 ];
    981     float stdev = myParams->data.F32[ 1 ];
     930float p_psGaussianDeriv(const psVector * restrict myData, const psVector * restrict myParams, int whichParam)
     931{
     932    float x = myData->data.F32[0];
     933    float mean = myParams->data.F32[0];
     934    float stdev = myParams->data.F32[1];
    982935    float tmp = 0.0;
    983936
    984     if ( whichParam == 0 ) {
     937    if (whichParam == 0) {
    985938        // Return the derivative w.r.t. the mean.
    986         tmp = ( x - mean ) * p_psGaussian( myData, myParams );
    987         tmp /= ( stdev * stdev );
    988     } else
    989         if ( whichParam == 1 ) {
    990             // Return the derivative w.r.t. the stdev.
    991             tmp = ( x - mean ) * ( x - mean ) * p_psGaussian( myData, myParams );
    992             tmp /= ( stdev * stdev * stdev );
    993         }
    994     printf( "p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp );
    995 
    996     return ( tmp );
    997 }
    998 
     939        tmp = (x - mean) * p_psGaussian(myData, myParams);
     940        tmp /= (stdev * stdev);
     941    } else if (whichParam == 1) {
     942        // Return the derivative w.r.t. the stdev.
     943        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
     944        tmp /= (stdev * stdev * stdev);
     945    }
     946    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
     947
     948    return (tmp);
     949}
    999950
    1000951/*****************************************************************************
    … …  
    1003954specified data point.
    1004955 *****************************************************************************/
    1005 float p_psQuadratic( const psVector *restrict myParams,
    1006                      const psVector *restrict myCoords )
    1007 {
    1008     float x = myCoords->data.F32[ 0 ];
    1009     float A = myParams->data.F32[ 0 ];
    1010     float B = myParams->data.F32[ 1 ];
    1011     float C = myParams->data.F32[ 2 ];
     956float p_psQuadratic(const psVector * restrict myParams, const psVector * restrict myCoords)
     957{
     958    float x = myCoords->data.F32[0];
     959    float A = myParams->data.F32[0];
     960    float B = myParams->data.F32[1];
     961    float C = myParams->data.F32[2];
    1012962    float tmp = 0.0;
    1013963
    1014     tmp = ( A * x * x ) + ( B * x ) + C;
    1015     return ( tmp );
     964    tmp = (A * x * x) + (B * x) + C;
     965    return (tmp);
    1016966}
    1017967
    … …  
    1020970calculates the specified partial derivative of the above quadratic function.
    1021971 *****************************************************************************/
    1022 float p_psQuadraticDeriv( const psVector *restrict myParams,
    1023                           const psVector *restrict myCoords,
    1024                           int whichParamDeriv )
    1025 {
    1026     float x = myCoords->data.F32[ 0 ];
     972float p_psQuadraticDeriv(const psVector * restrict myParams,
     973                         const psVector * restrict myCoords, int whichParamDeriv)
     974{
     975    float x = myCoords->data.F32[0];
    1027976    float tmp = 0.0;
    1028977
    1029     if ( whichParamDeriv == 0 ) {
     978    if (whichParamDeriv == 0) {
    1030979        tmp = x * x;
    1031     } else
    1032         if ( whichParamDeriv == 1 ) {
    1033             tmp = x;
    1034         } else
    1035             if ( whichParamDeriv == 2 ) {
    1036                 tmp = 1.0;
    1037             }
    1038 
    1039     return ( tmp );
     980    } else if (whichParamDeriv == 1) {
     981        tmp = x;
     982    } else if (whichParamDeriv == 2) {
     983        tmp = 1.0;
     984    }
     985
     986    return (tmp);
    1040987}
    1041988
    … …  
    1049996decreasing within that range.
    1050997 *****************************************************************************/
    1051 float p_ps1DPolyMedian( psPolynomial1D *myPoly,
    1052                         float rangeLow,
    1053                         float rangeHigh,
    1054                         float getThisValue )
     998float p_ps1DPolyMedian(psPolynomial1D * myPoly, float rangeLow, float rangeHigh, float getThisValue)
    1055999{
    10561000    int numIterations = 0;
    … …  
    10591003    float f = 0.0;
    10601004
    1061     //  printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
    1062 
    1063     while ( numIterations < MAX_ITERATIONS ) {
    1064         midpoint = ( rangeHigh + rangeLow ) / 2.0;
    1065         if ( fabs( midpoint - oldMidpoint ) <= FLT_EPSILON ) {
    1066             return ( midpoint );
     1005    // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
     1006
     1007    while (numIterations < MAX_ITERATIONS) {
     1008        midpoint = (rangeHigh + rangeLow) / 2.0;
     1009        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
     1010            return (midpoint);
    10671011        }
    10681012        oldMidpoint = midpoint;
    10691013
    1070         f = psPolynomial1DEval( midpoint, myPoly );
    1071         //        printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
    1072         if ( fabs( f - getThisValue ) <= FLT_EPSILON ) {
    1073             return ( midpoint );
    1074         }
    1075 
    1076         if ( f > getThisValue ) {
     1014        f = psPolynomial1DEval(midpoint, myPoly);
     1015        // printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
     1016        if (fabs(f - getThisValue) <= FLT_EPSILON) {
     1017            return (midpoint);
     1018        }
     1019
     1020        if (f > getThisValue) {
    10771021            rangeHigh = midpoint;
    10781022        } else {
    … …  
    10811025        numIterations++;
    10821026    }
    1083     return ( midpoint );
     1027    return (midpoint);
    10841028}
    10851029
    … …  
    10921036XXX: This function is currently not being used.
    10931037 *****************************************************************************/
    1094 float p_psFitQuadratic( psHistogram *histogram,
    1095                         psVector *cumulativeSums,
    1096                         int binNum,
    1097                         float fitFloat )
    1098 {
    1099     psVector * x = psVectorAlloc( 3, PS_TYPE_F64 );
    1100     psVector *y = psVectorAlloc( 3, PS_TYPE_F64 );
    1101     psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 );
    1102     psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 );
    1103 
    1104     if ( ( binNum > 0 ) &&
    1105             ( binNum < ( histogram->nums->n + 1 ) ) ) {
    1106         x->data.F64[ 0 ] = ( double ) 0.5 *
    1107                            ( histogram->bounds->data.F32[ binNum - 1 ] +
    1108                              histogram->bounds->data.F32[ binNum ] );
    1109         x->data.F64[ 1 ] = ( double ) 0.5 *
    1110                            ( histogram->bounds->data.F32[ binNum ] +
    1111                              histogram->bounds->data.F32[ binNum + 1 ] );
    1112         x->data.F64[ 2 ] = ( double ) 0.5 *
    1113                            ( histogram->bounds->data.F32[ binNum + 1 ] +
    1114                              histogram->bounds->data.F32[ binNum + 2 ] );
    1115 
    1116         y->data.F64[ 0 ] = cumulativeSums->data.F32[ binNum - 1 ];
    1117         y->data.F64[ 1 ] = cumulativeSums->data.F32[ binNum ];
    1118         y->data.F64[ 2 ] = cumulativeSums->data.F32[ binNum + 1 ];
    1119 
    1120         if ( !( ( y->data.F64[ 0 ] <= fitFloat ) &&
    1121                 ( fitFloat <= y->data.F64[ 2 ] ) ) ) {
    1122             psAbort( __func__, "p_psVectorRobustStats(0): midpoint not within y-range\n" );
    1123         }
    1124 
    1125         yErr->data.F64[ 0 ] = 1.0;
    1126         yErr->data.F64[ 1 ] = 1.0;
    1127         yErr->data.F64[ 2 ] = 1.0;
    1128 
    1129         myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
    1130         return ( p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ], x->data.F64[ 2 ],
    1131                                    fitFloat ) );
     1038float p_psFitQuadratic(psHistogram * histogram, psVector * cumulativeSums, int binNum, float fitFloat)
     1039{
     1040    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1041    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1042    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1043    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
     1044
     1045    if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) {
     1046        x->data.F64[0] = (double)0.5 *
     1047                         (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]);
     1048        x->data.F64[1] = (double)0.5 *
     1049                         (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]);
     1050        x->data.F64[2] = (double)0.5 *
     1051                         (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]);
     1052
     1053        y->data.F64[0] = cumulativeSums->data.F32[binNum - 1];
     1054        y->data.F64[1] = cumulativeSums->data.F32[binNum];
     1055        y->data.F64[2] = cumulativeSums->data.F32[binNum + 1];
     1056
     1057        if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) {
     1058            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
     1059        }
     1060
     1061        yErr->data.F64[0] = 1.0;
     1062        yErr->data.F64[1] = 1.0;
     1063        yErr->data.F64[2] = 1.0;
     1064
     1065        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
     1066        return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat));
    11321067    } else {
    1133         return ( 0.5 * ( histogram->bounds->data.F32[ binNum + 1 ] +
    1134                          histogram->bounds->data.F32[ binNum ] ) );
    1135     }
    1136 
    1137     psFree( x );
    1138     psFree( y );
    1139     psFree( yErr );
    1140     psFree( myPoly );
    1141     return ( 0.0 );
     1068        return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum]));
     1069    }
     1070
     1071    psFree(x);
     1072    psFree(y);
     1073    psFree(yErr);
     1074    psFree(myPoly);
     1075    return (0.0);
    11421076}
    11431077
    … …  
    11661100    NULL
    11671101*****************************************************************************/
    1168 void p_psVectorRobustStats( const psVector *restrict myVector,
    1169                             const psVector *restrict maskVector,
    1170                             unsigned int maskVal,
    1171                             psStats *stats )
    1172 {
    1173     psHistogram * robustHistogram = NULL;
     1102void p_psVectorRobustStats(const psVector * restrict myVector,
     1103                           const psVector * restrict maskVector, unsigned int maskVal, psStats * stats)
     1104{
     1105    psHistogram *robustHistogram = NULL;
    11741106    psVector *robustHistogramVector = NULL;
    1175     float binSize = 0.0;                        // Size of the histogram bins
    1176     int LQBinNum = -1;                          // Bin num for lower quartile
    1177     int UQBinNum = -1;                          // Bin num for upper quartile
    1178     int i = 0;                                  // Loop index variable
     1107    float binSize = 0.0;        // Size of the histogram bins
     1108    int LQBinNum = -1;          // Bin num for lower quartile
     1109    int UQBinNum = -1;          // Bin num for upper quartile
     1110    int i = 0;                  // Loop index variable
    11791111    int maxBinNum = 0;
    11801112    float maxBinCount = 0.0;
    11811113    float dL = 0.0;
    11821114    int numBins = 0;
    1183     psStats *tmpStats = psStatsAlloc( PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN );
    1184     //    psImage *domain;
    1185     //    psVector *errors;
    1186     //    psVector *data;
    1187     //    psVector *initialGuess;
    1188     //    psVector *theParams;
    1189     //    float chiSq=0.0;
    1190     //    float max = -HUGE;
    1191     //    float max_pos;
     1115    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
     1116
     1117    // psImage *domain;
     1118    // psVector *errors;
     1119    // psVector *data;
     1120    // psVector *initialGuess;
     1121    // psVector *theParams;
     1122    // float chiSq=0.0;
     1123    // float max = -HUGE;
     1124    // float max_pos;
    11921125    float myMean = 0.0;
    11931126    float myStdev = 0.0;
    … …  
    11961129    float sumSquares = 0.0;
    11971130    float sumDiffs = 0.0;
    1198     psVector *x = psVectorAlloc( 3, PS_TYPE_F64 );
    1199     psVector *y = psVectorAlloc( 3, PS_TYPE_F64 );
    1200     psVector *yErr = psVectorAlloc( 3, PS_TYPE_F64 );
    1201     psPolynomial1D *myPoly = psPolynomial1DAlloc( 2 );
     1131    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1132    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1133    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1134    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
    12021135    psVector *cumulativeRobustSumsFullRange = NULL;
    12031136    psVector *cumulativeRobustSumsDlRange = NULL;
    … …  
    12101143    // by computing the clipped standard deviation of the vector, and dividing
    12111144    // that by 10.0;
    1212     p_psVectorClippedStats( myVector, maskVector, maskVal, tmpStats );
     1145    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
    12131146    binSize = tmpStats->clippedStdev / 10.0f;
    12141147
    12151148    // If stats->clippedStdev == 0.0, then all data elements have the same
    12161149    // value.  Therefore, we can set the appropiate results and return.
    1217     if ( fabs( binSize ) <= FLT_EPSILON ) {
    1218         if ( stats->options & PS_STAT_ROBUST_MEAN ) {
     1150    if (fabs(binSize) <= FLT_EPSILON) {
     1151        if (stats->options & PS_STAT_ROBUST_MEAN) {
    12191152            stats->robustMean = stats->clippedMean;
    12201153        }
    1221         if ( stats->options & PS_STAT_ROBUST_MEDIAN ) {
     1154        if (stats->options & PS_STAT_ROBUST_MEDIAN) {
    12221155            stats->robustMedian = stats->clippedMean;
    12231156        }
    1224         if ( stats->options & PS_STAT_ROBUST_MODE ) {
     1157        if (stats->options & PS_STAT_ROBUST_MODE) {
    12251158            stats->robustMode = stats->clippedMean;
    12261159        }
    1227         if ( stats->options & PS_STAT_ROBUST_STDEV ) {
     1160        if (stats->options & PS_STAT_ROBUST_STDEV) {
    12281161            stats->robustStdev = 0.0;
    12291162        }
    1230         if ( stats->options & PS_STAT_ROBUST_QUARTILE ) {
     1163        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
    12311164            stats->robustUQ = stats->clippedMean;
    12321165            stats->robustLQ = stats->clippedMean;
    … …  
    12351168        stats->robustNfit = 0.0;
    12361169        stats->robustN50 = 0.0;
    1237         psFree( tmpStats );
    1238         return ;
    1239     }
    1240 
     1170        psFree(tmpStats);
     1171        return;
     1172    }
    12411173    // Determine minimum and maximum values in the data vector.
    1242     if ( isnan( stats->min ) ) {
    1243         p_psVectorMin( myVector, maskVector, maskVal, stats );
    1244     }
    1245     if ( isnan( stats->max ) ) {
    1246         p_psVectorMax( myVector, maskVector, maskVal, stats );
    1247     }
    1248 
     1174    if (isnan(stats->min)) {
     1175        p_psVectorMin(myVector, maskVector, maskVal, stats);
     1176    }
     1177    if (isnan(stats->max)) {
     1178        p_psVectorMax(myVector, maskVector, maskVal, stats);
     1179    }
    12491180    // Create the histogram structure.  NOTE: we can not specify the bin size
    12501181    // precisely since the argument to psHistogramAlloc() is the number of
    12511182    // bins, not the binSize.  Also, if we get here, we know that
    12521183    // binSize != 0.0.
    1253     numBins = ( int ) ( ( stats->max - stats->min ) / binSize );
    1254     robustHistogram = psHistogramAlloc( stats->min,
    1255                                         stats->max,
    1256                                         numBins );
     1184    numBins = (int)((stats->max - stats->min) / binSize);
     1185    robustHistogram = psHistogramAlloc(stats->min, stats->max, numBins);
    12571186
    12581187    // Populate the histogram array.
    1259     psVectorHistogram( robustHistogram, myVector, maskVector, maskVal );
     1188    psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
    12601189
    12611190    // Smooth the histogram.
    12621191    // XXX: is that the right stdev?
    1263     robustHistogramVector = p_psVectorsmoothHistGaussian( robustHistogram,
    1264                             tmpStats->clippedStdev / 4.0f );
     1192    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, tmpStats->clippedStdev / 4.0f);
    12651193
    12661194    // The following was necessary to fit a gaussian to the data, since
    … …  
    12741202    // index position i is equal to the sum of bins 0:i.  This will be used
    12751203    // now and later in determining the lower/upper quartiles.
    1276     cumulativeRobustSumsFullRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 );
    1277     cumulativeRobustSumsFullRange->data.F32[ 0 ] = robustHistogramVector->data.F32[ 0 ];
    1278     for ( i = 1;i < robustHistogramVector->n;i++ ) {
    1279         cumulativeRobustSumsFullRange->data.F32[ i ] =
    1280             cumulativeRobustSumsFullRange->data.F32[ i - 1 ] +
    1281             robustHistogramVector->data.F32[ i ];
    1282     }
    1283     sumRobust = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ];
     1204    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1205    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
     1206    for (i = 1; i < robustHistogramVector->n; i++) {
     1207        cumulativeRobustSumsFullRange->data.F32[i] =
     1208            cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
     1209    }
     1210    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
    12841211
    12851212    // Determine the bin number containing the lower quartile point.
    12861213    LQBinNum = -1;
    1287     for ( i = 0;i < cumulativeRobustSumsFullRange->n;i++ ) {
    1288         if ( cumulativeRobustSumsFullRange->data.F32[ i ] >= ( sumRobust / 4.0 ) ) {
     1214    for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) {
     1215        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) {
    12891216            LQBinNum = i;
    12901217            break;
    … …  
    12941221    // Determine the bin number containing the upper quartile point.
    12951222    UQBinNum = -1;
    1296     for ( i = cumulativeRobustSumsFullRange->n - 1;i >= 0;i-- ) {
    1297         if ( cumulativeRobustSumsFullRange->data.F32[ i ] <= ( 3.0 * sumRobust / 4.0 ) ) {
     1223    for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) {
     1224        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) {
    12981225            UQBinNum = i;
    12991226            break;
    … …  
    13011228    }
    13021229
    1303     if ( ( LQBinNum == -1 ) ||
    1304             ( UQBinNum == -1 ) ) {
    1305         psAbort( __func__, "Could not determine the robust lower/upper quartiles." );
    1306     }
     1230    if ((LQBinNum == -1) || (UQBinNum == -1)) {
     1231        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
     1232    }
     1233
    13071234    /**************************************************************************
    13081235    Determine the mode in the range LQ:UQ.
    … …  
    13101237    // Determine the bin with the peak value in the range LQ to UQ.
    13111238    maxBinNum = LQBinNum;
    1312     maxBinCount = robustHistogramVector->data.F32[ LQBinNum ];
    1313     sumN50 = ( float ) robustHistogram->nums->data.S32[ LQBinNum ];
    1314     for ( i = LQBinNum + 1;i <= UQBinNum;i++ ) {
    1315         if ( robustHistogramVector->data.F32[ i ] > maxBinCount ) {
     1239    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
     1240    sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum];
     1241    for (i = LQBinNum + 1; i <= UQBinNum; i++) {
     1242        if (robustHistogramVector->data.F32[i] > maxBinCount) {
    13161243            maxBinNum = i;
    1317             maxBinCount = robustHistogramVector->data.F32[ i ];
    1318         }
    1319         sumN50 += ( float ) robustHistogram->nums->data.S32[ i ];
     1244            maxBinCount = robustHistogramVector->data.F32[i];
     1245        }
     1246        sumN50 += (float)robustHistogram->nums->data.S32[i];
    13201247    }
    13211248
    13221249    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
    1323     dL = ( UQBinNum - LQBinNum ) / 4;
    1324 
    1325     printf( "(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum );
     1250    dL = (UQBinNum - LQBinNum) / 4;
     1251
     1252    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
    13261253
    13271254    /**************************************************************************
    13281255    Determine the mean/stdev for the bins in the range mode-dL to mode+dL
    13291256    **************************************************************************/
    1330     cumulativeRobustSumsDlRange = psVectorAlloc( robustHistogramVector->n, PS_TYPE_F32 );
    1331     for ( i = 0;i < robustHistogramVector->n;i++ ) {
    1332         cumulativeRobustSumsDlRange->data.F32[ i ] = 0.0;
     1257    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1258    for (i = 0; i < robustHistogramVector->n; i++) {
     1259        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
    13331260    }
    13341261    sumNfit = 0.0;
    13351262    cumulativeMedian = 0.0;
    1336     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
    1337         if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
    1338             cumulativeRobustSumsDlRange->data.F32[ i ] =
    1339                 cumulativeRobustSumsDlRange->data.F32[ i - 1 ] +
    1340                 robustHistogramVector->data.F32[ i ];
    1341             cumulativeMedian += robustHistogramVector->data.F32[ i ];
    1342             sumNfit += ( float ) robustHistogram->nums->data.S32[ i ];
     1263    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
     1264        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1265            cumulativeRobustSumsDlRange->data.F32[i] =
     1266                cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
     1267            cumulativeMedian += robustHistogramVector->data.F32[i];
     1268            sumNfit += (float)robustHistogram->nums->data.S32[i];
    13431269        }
    13441270    }
    … …  
    13481274    // that bin (this is a non-exact approximation).
    13491275    myMean = 0.0;
    1350     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
    1351         if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
    1352             myMean += ( robustHistogramVector->data.F32[ i ] ) * 0.5 *
    1353                       ( robustHistogram->bounds->data.F32[ i + 1 ] +
    1354                         robustHistogram->bounds->data.F32[ i ] );
    1355             countFloat += robustHistogramVector->data.F32[ i ];
     1276    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
     1277        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1278            myMean += (robustHistogramVector->data.F32[i]) * 0.5 *
     1279                      (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]);
     1280            countFloat += robustHistogramVector->data.F32[i];
    13561281        }
    13571282    }
    … …  
    13611286    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
    13621287    // that bin.
    1363     for ( i = maxBinNum - dL;i <= maxBinNum + dL;i++ ) {
    1364         if ( ( 0 <= i ) && ( i < robustHistogramVector->n ) ) {
    1365             diff = ( 0.5 * ( robustHistogram->bounds->data.F32[ i + 1 ] +
    1366                              robustHistogram->bounds->data.F32[ i ] ) ) - myMean;
    1367             sumSquares += diff * diff * robustHistogramVector->data.F32[ i ];
    1368             sumDiffs += diff * robustHistogramVector->data.F32[ i ];
    1369         }
    1370     }
    1371     myStdev = sqrt( ( sumSquares - ( sumDiffs * sumDiffs / countFloat ) ) / ( countFloat - 1 ) );
     1288    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
     1289        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1290            diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] +
     1291                           robustHistogram->bounds->data.F32[i])) - myMean;
     1292            sumSquares += diff * diff * robustHistogramVector->data.F32[i];
     1293            sumDiffs += diff * robustHistogramVector->data.F32[i];
     1294        }
     1295    }
     1296    myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    13721297
    13731298    /**************************************************************************
    13741299    Set the appropriate members in the output stats struct.
    13751300    **************************************************************************/
    1376     if ( stats->options & PS_STAT_ROBUST_MEAN ) {
     1301    if (stats->options & PS_STAT_ROBUST_MEAN) {
    13771302        stats->robustMean = myMean;
    13781303    }
    13791304
    1380     if ( stats->options & PS_STAT_ROBUST_MODE ) {
     1305    if (stats->options & PS_STAT_ROBUST_MODE) {
    13811306        stats->robustMode = 0.5 *
    1382                             ( robustHistogram->bounds->data.F32[ maxBinNum ] +
    1383                               robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
    1384     }
    1385 
    1386     if ( stats->options & PS_STAT_ROBUST_STDEV ) {
     1307                            (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]);
     1308    }
     1309
     1310    if (stats->options & PS_STAT_ROBUST_STDEV) {
    13871311        stats->robustStdev = myStdev;
    13881312    }
    1389 
    13901313    // To determine the median (and later, the lower/upper quartile), we fit
    13911314    // a quadratic to the three bins surrounding the bin containing the median.
    … …  
    13941317    // this bin.  We then solve the quadratic for
    13951318
    1396     if ( stats->options & PS_STAT_ROBUST_MEDIAN ) {
    1397         if ( ( maxBinNum > 0 ) && ( maxBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
    1398             x->data.F64[ 0 ] = ( double ) 0.5 *
    1399                                ( robustHistogram->bounds->data.F32[ maxBinNum - 1 ] +
    1400                                  robustHistogram->bounds->data.F32[ maxBinNum ] );
    1401             x->data.F64[ 1 ] = ( double ) 0.5 *
    1402                                ( robustHistogram->bounds->data.F32[ maxBinNum ] +
    1403                                  robustHistogram->bounds->data.F32[ maxBinNum + 1 ] );
    1404             x->data.F64[ 2 ] = ( double ) 0.5 *
    1405                                ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
    1406                                  robustHistogram->bounds->data.F32[ maxBinNum + 2 ] );
    1407 
    1408             y->data.F64[ 0 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum - 1 ];
    1409             y->data.F64[ 1 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum ];
    1410             y->data.F64[ 2 ] = cumulativeRobustSumsDlRange->data.F32[ maxBinNum + 1 ];
     1319    if (stats->options & PS_STAT_ROBUST_MEDIAN) {
     1320        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) {
     1321            x->data.F64[0] = (double)0.5 *
     1322                             (robustHistogram->bounds->data.F32[maxBinNum - 1] +
     1323                              robustHistogram->bounds->data.F32[maxBinNum]);
     1324            x->data.F64[1] = (double)0.5 *
     1325                             (robustHistogram->bounds->data.F32[maxBinNum] +
     1326                              robustHistogram->bounds->data.F32[maxBinNum + 1]);
     1327            x->data.F64[2] = (double)0.5 *
     1328                             (robustHistogram->bounds->data.F32[maxBinNum + 1] +
     1329                              robustHistogram->bounds->data.F32[maxBinNum + 2]);
     1330
     1331            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1];
     1332            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
     1333            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1];
    14111334
    14121335            // Ensure that cumulativeMedian/2 is actually within the range of the bins
    14131336            // we are using.
    14141337            cumulativeMedian *= 0.5;
    1415             if ( !( ( y->data.F64[ 0 ] <= cumulativeMedian ) &&
    1416                     ( cumulativeMedian <= y->data.F64[ 2 ] ) ) ) {
    1417                 printf( "((%f), %f, %f)\n", cumulativeMedian, y->data.F64[ 0 ], y->data.F64[ 2 ] );
    1418                 psAbort( __func__, "p_psVectorRobustStats(1): midpoint not within y-range\n" );
     1338            if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) {
     1339                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
     1340                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
    14191341            }
    14201342            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
    14211343            // may have to set this meaningfully later.
    1422             yErr->data.F64[ 0 ] = 1.0;
    1423             yErr->data.F64[ 1 ] = 1.0;
    1424             yErr->data.F64[ 2 ] = 1.0;
     1344            yErr->data.F64[0] = 1.0;
     1345            yErr->data.F64[1] = 1.0;
     1346            yErr->data.F64[2] = 1.0;
    14251347
    14261348            // Determine the coefficients of the polynomial.
    1427             myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
     1349            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
    14281350            // Call p_ps1DPolyMedian(), which does a binary search on the
    14291351            // polynomial, looking for the value x such that
    14301352            // f(x) = cumulativeMedian.
    1431             stats->robustMedian = p_ps1DPolyMedian( myPoly, x->data.F64[ 0 ],
    1432                                                     x->data.F64[ 2 ], cumulativeMedian );
     1353            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian);
    14331354        } else {
    14341355            // If the mode is the first/last histogram bin, then simply use
    14351356            // the midpoint of that bin.
    1436             stats->robustMedian = 0.5 * ( robustHistogram->bounds->data.F32[ maxBinNum + 1 ] +
    1437                                           robustHistogram->bounds->data.F32[ maxBinNum ] );
    1438         }
    1439     }
    1440 
     1357            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] +
     1358                                         robustHistogram->bounds->data.F32[maxBinNum]);
     1359        }
     1360    }
    14411361    // The lower/upper quartile calculations are very similar to the median
    14421362    // calculations.  We fit a quadratic to the array containing the
    … …  
    14441364    // f(x) equals the lower/upper quartile exactly.
    14451365    //
    1446     if ( stats->options & PS_STAT_ROBUST_QUARTILE ) {
    1447         countFloat = cumulativeRobustSumsFullRange->data.F32[ robustHistogramVector->n - 1 ];
    1448 
    1449         if ( ( LQBinNum > 0 ) && ( LQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
    1450             x->data.F64[ 0 ] = ( double ) 0.5 *
    1451                                ( robustHistogram->bounds->data.F32[ LQBinNum - 1 ] +
    1452                                  robustHistogram->bounds->data.F32[ LQBinNum ] );
    1453             x->data.F64[ 1 ] = ( double ) 0.5 *
    1454                                ( robustHistogram->bounds->data.F32[ LQBinNum ] +
    1455                                  robustHistogram->bounds->data.F32[ LQBinNum + 1 ] );
    1456             x->data.F64[ 2 ] = ( double ) 0.5 *
    1457                                ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
    1458                                  robustHistogram->bounds->data.F32[ LQBinNum + 2 ] );
    1459 
    1460             y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum - 1 ];
    1461             y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum ];
    1462             y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ LQBinNum + 1 ];
    1463 
    1464             if ( !( ( y->data.F64[ 0 ] <= ( countFloat / 4.0 ) ) &&
    1465                     ( ( countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
    1466                 psAbort( __func__, "p_psVectorRobustStats(2): midpoint not within y-range\n" );
    1467             }
    1468 
    1469             yErr->data.F64[ 0 ] = 1.0;
    1470             yErr->data.F64[ 1 ] = 1.0;
    1471             yErr->data.F64[ 2 ] = 1.0;
    1472 
    1473             myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
    1474             stats->robustLQ = p_ps1DPolyMedian( myPoly,
    1475                                                 x->data.F64[ 0 ],
    1476                                                 x->data.F64[ 2 ],
    1477                                                 countFloat / 4.0 );
     1366    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
     1367        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
     1368
     1369        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) {
     1370            x->data.F64[0] = (double)0.5 *
     1371                             (robustHistogram->bounds->data.F32[LQBinNum - 1] +
     1372                              robustHistogram->bounds->data.F32[LQBinNum]);
     1373            x->data.F64[1] = (double)0.5 *
     1374                             (robustHistogram->bounds->data.F32[LQBinNum] +
     1375                              robustHistogram->bounds->data.F32[LQBinNum + 1]);
     1376            x->data.F64[2] = (double)0.5 *
     1377                             (robustHistogram->bounds->data.F32[LQBinNum + 1] +
     1378                              robustHistogram->bounds->data.F32[LQBinNum + 2]);
     1379
     1380            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1];
     1381            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
     1382            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1];
     1383
     1384            if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) {
     1385                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
     1386            }
     1387
     1388            yErr->data.F64[0] = 1.0;
     1389            yErr->data.F64[1] = 1.0;
     1390            yErr->data.F64[2] = 1.0;
     1391
     1392            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
     1393            stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0);
    14781394
    14791395        } else {
    14801396            // If the LQ is the first/last histogram bin, then simply use
    14811397            // the midpoint of that bin.
    1482             stats->robustLQ = 0.5 * ( robustHistogram->bounds->data.F32[ LQBinNum + 1 ] +
    1483                                       robustHistogram->bounds->data.F32[ LQBinNum ] );
    1484         }
    1485 
    1486         if ( ( UQBinNum > 0 ) && ( UQBinNum < ( robustHistogram->nums->n - 1 ) ) ) {
    1487             x->data.F64[ 0 ] = ( double ) 0.5 *
    1488                                ( robustHistogram->bounds->data.F32[ UQBinNum - 1 ] +
    1489                                  robustHistogram->bounds->data.F32[ UQBinNum ] );
    1490             x->data.F64[ 1 ] = ( double ) 0.5 *
    1491                                ( robustHistogram->bounds->data.F32[ UQBinNum ] +
    1492                                  robustHistogram->bounds->data.F32[ UQBinNum + 1 ] );
    1493             x->data.F64[ 2 ] = ( double ) 0.5 *
    1494                                ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
    1495                                  robustHistogram->bounds->data.F32[ UQBinNum + 2 ] );
    1496 
    1497             y->data.F64[ 0 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum - 1 ];
    1498             y->data.F64[ 1 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum ];
    1499             y->data.F64[ 2 ] = cumulativeRobustSumsFullRange->data.F32[ UQBinNum + 1 ];
    1500 
    1501             if ( !( ( y->data.F64[ 0 ] <= ( 3.0 * countFloat / 4.0 ) ) &&
    1502                     ( ( 3.0 * countFloat / 4.0 ) <= y->data.F64[ 2 ] ) ) ) {
    1503                 psAbort( __func__, "p_psVectorRobustStats(3): midpoint not within y-range\n" );
    1504             }
    1505 
    1506             yErr->data.F64[ 0 ] = 1.0;
    1507             yErr->data.F64[ 1 ] = 1.0;
    1508             yErr->data.F64[ 2 ] = 1.0;
    1509 
    1510             myPoly = psVectorFitPolynomial1D( myPoly, x, y, yErr );
    1511             stats->robustUQ = p_ps1DPolyMedian( myPoly,
    1512                                                 x->data.F64[ 0 ],
    1513                                                 x->data.F64[ 2 ],
    1514                                                 3.0 * countFloat / 4.0 );
     1398            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] +
     1399                                     robustHistogram->bounds->data.F32[LQBinNum]);
     1400        }
     1401
     1402        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) {
     1403            x->data.F64[0] = (double)0.5 *
     1404                             (robustHistogram->bounds->data.F32[UQBinNum - 1] +
     1405                              robustHistogram->bounds->data.F32[UQBinNum]);
     1406            x->data.F64[1] = (double)0.5 *
     1407                             (robustHistogram->bounds->data.F32[UQBinNum] +
     1408                              robustHistogram->bounds->data.F32[UQBinNum + 1]);
     1409            x->data.F64[2] = (double)0.5 *
     1410                             (robustHistogram->bounds->data.F32[UQBinNum + 1] +
     1411                              robustHistogram->bounds->data.F32[UQBinNum + 2]);
     1412
     1413            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1];
     1414            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
     1415            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1];
     1416
     1417            if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) &&
     1418                    ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) {
     1419                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
     1420            }
     1421
     1422            yErr->data.F64[0] = 1.0;
     1423            yErr->data.F64[1] = 1.0;
     1424            yErr->data.F64[2] = 1.0;
     1425
     1426            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
     1427            stats->robustUQ = p_ps1DPolyMedian(myPoly,
     1428                                               x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0);
    15151429        } else {
    15161430            // If the UQ is the first/last histogram bin, then simply use
    15171431            // the midpoint of that bin.
    1518             stats->robustUQ = 0.5 * ( robustHistogram->bounds->data.F32[ UQBinNum + 1 ] +
    1519                                       robustHistogram->bounds->data.F32[ UQBinNum ] );
     1432            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] +
     1433                                     robustHistogram->bounds->data.F32[UQBinNum]);
    15201434        }
    15211435    }
    … …  
    15231437    stats->robustN50 = sumN50;
    15241438
    1525     psFree( x );
    1526     psFree( y );
    1527     psFree( yErr );
    1528     psFree( tmpStats );
    1529     psFree( robustHistogram );
    1530     psFree( myPoly );
    1531     psFree( cumulativeRobustSumsFullRange );
    1532     psFree( cumulativeRobustSumsDlRange );
    1533 }
    1534 
    1535 
     1439    psFree(x);
     1440    psFree(y);
     1441    psFree(yErr);
     1442    psFree(tmpStats);
     1443    psFree(robustHistogram);
     1444    psFree(myPoly);
     1445    psFree(cumulativeRobustSumsFullRange);
     1446    psFree(cumulativeRobustSumsDlRange);
     1447}
    15361448
    15371449/*
    … …  
    15821494*/
    15831495
    1584 
    15851496/*****************************************************************************/
     1497
    15861498/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     1499
    15871500/*****************************************************************************/
    15881501
    1589 static void histogramFree( psHistogram *myHist );
     1502static void histogramFree(psHistogram * myHist);
    15901503
    15911504/******************************************************************************
    15921505    psStatsAlloc(): This routine must create a new psStats data structure.
    15931506 *****************************************************************************/
    1594 psStats *psStatsAlloc( psStatsOptions options )
    1595 {
    1596     psStats * newStruct = NULL;
    1597 
    1598     newStruct = ( psStats * ) psAlloc( sizeof( psStats ) );
     1507psStats *psStatsAlloc(psStatsOptions options)
     1508{
     1509    psStats *newStruct = NULL;
     1510
     1511    newStruct = (psStats *) psAlloc(sizeof(psStats));
    15991512    newStruct->sampleMean = NAN;
    16001513    newStruct->sampleMedian = NAN;
    … …  
    16201533    newStruct->options = options;
    16211534
    1622     return ( newStruct );
     1535    return (newStruct);
    16231536}
    16241537
    … …  
    16351548    The histogram structure
    16361549 *****************************************************************************/
    1637 psHistogram *psHistogramAlloc( float lower,
    1638                                float upper,
    1639                                int n )
    1640 {
    1641     int i = 0;                          // Loop index variable
     1550psHistogram *psHistogramAlloc(float lower, float upper, int n)
     1551{
     1552    int i = 0;                  // Loop index variable
    16421553    psHistogram *newHist = NULL;        // The new histogram structure
    1643     float binSize = 0.0;                // The histogram bin size
     1554    float binSize = 0.0;        // The histogram bin size
    16441555
    16451556    // NOTE: Verify that this is the correct action.
    1646     if ( n == 0 ) {
    1647         return ( NULL );
    1648     }
    1649     if ( n < 0 ) {
    1650         psAbort( __func__, "psHistogramAlloc() called with bin size %d.\n", n );
    1651     }
    1652 
     1557    if (n == 0) {
     1558        return (NULL);
     1559    }
     1560    if (n < 0) {
     1561        psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n);
     1562    }
    16531563    // NOTE: Verify that this is the correct action.
    1654     if ( lower > upper ) {
    1655         return ( NULL );
    1656     }
    1657 
     1564    if (lower > upper) {
     1565        return (NULL);
     1566    }
    16581567    // Allocate memory for the new histogram structure.  If there are N
    16591568    // bins, then there are N+1 bounds to those bins.
    1660     newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) );
    1661     p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree );
    1662     newHist->bounds = psVectorAlloc( n + 1, PS_TYPE_F32 );
     1569    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     1570    p_psMemSetDeallocator(newHist, (psFreeFcn) histogramFree);
     1571    newHist->bounds = psVectorAlloc(n + 1, PS_TYPE_F32);
    16631572    newHist->bounds->n = newHist->bounds->nalloc;
    16641573
    16651574    // Calculate the bounds for each bin.
    1666     binSize = ( upper - lower ) / ( float ) n;
    1667     // NOTE: Is the following necessary?  It prevents the max data point
     1575    binSize = (upper - lower) / (float)n;
     1576    // NOTE: Is the following necessary? It prevents the max data point
    16681577    // from being in a non-existant bin.
    16691578    binSize += FLT_EPSILON;
    1670     for ( i = 0;i < n + 1;i++ ) {
    1671         newHist->bounds->data.F32[ i ] = lower + ( binSize * ( float ) i );
     1579    for (i = 0; i < n + 1; i++) {
     1580        newHist->bounds->data.F32[i] = lower + (binSize * (float)i);
    16721581    }
    16731582
    16741583    // Allocate the bins, and initialize them to zero.
    1675     newHist->nums = psVectorAlloc( n, PS_TYPE_U32 );
     1584    newHist->nums = psVectorAlloc(n, PS_TYPE_U32);
    16761585    newHist->nums->n = newHist->nums->nalloc;
    1677     for ( i = 0;i < newHist->nums->n;i++ ) {
    1678         newHist->nums->data.U32[ i ] = 0;
     1586    for (i = 0; i < newHist->nums->n; i++) {
     1587        newHist->nums->data.U32[i] = 0;
    16791588    }
    16801589
    … …  
    16841593    newHist->uniform = true;
    16851594
    1686     return ( newHist );
     1595    return (newHist);
    16871596}
    16881597
    … …  
    16961605    The histogram structure
    16971606 *****************************************************************************/
    1698 psHistogram *psHistogramAllocGeneric( const psVector *restrict bounds )
    1699 {
    1700     psHistogram * newHist = NULL;        // The new histogram structure
    1701     int i;                              // Loop index variable
     1607psHistogram *psHistogramAllocGeneric(const psVector * restrict bounds)
     1608{
     1609    psHistogram *newHist = NULL;        // The new histogram structure
     1610    int i;                      // Loop index variable
    17021611
    17031612    // NOTE: Verify that this is the correct action.
    1704     if ( bounds == NULL ) {
     1613    if (bounds == NULL) {
    17051614        // psAbort(__func__, "psHistogram requested with NULL bounds");
    1706         return ( NULL );
    1707     }
    1708 
     1615        return (NULL);
     1616    }
    17091617    // NOTE: Verify that this is the correct action.
    1710     if ( bounds->n <= 1 ) {
     1618    if (bounds->n <= 1) {
    17111619        // psAbort(__func__, "psHistogram requested with NULL bounds");
    1712         return ( NULL );
    1713     }
    1714 
    1715     if ( bounds->type.type != PS_TYPE_F32 ) {
     1620        return (NULL);
     1621    }
     1622
     1623    if (bounds->type.type != PS_TYPE_F32) {
    17161624        // psAbort(__func__, "psHistogram request a bound which is not type F32");
    1717         return ( NULL );
    1718     }
    1719 
     1625        return (NULL);
     1626    }
    17201627    // Allocate memory for the new histogram structure.
    1721     newHist = ( psHistogram * ) psAlloc( sizeof( psHistogram ) );
    1722     p_psMemSetDeallocator( newHist, ( psFreeFcn ) histogramFree );
    1723     newHist->bounds = psVectorAlloc( bounds->n, PS_TYPE_F32 );
     1628    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
     1629    p_psMemSetDeallocator(newHist, (psFreeFcn) histogramFree);
     1630    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
    17241631    newHist->bounds->n = newHist->bounds->nalloc;
    1725     for ( i = 0;i < bounds->n;i++ ) {
    1726         newHist->bounds->data.F32[ i ] = bounds->data.F32[ i ];
     1632    for (i = 0; i < bounds->n; i++) {
     1633        newHist->bounds->data.F32[i] = bounds->data.F32[i];
    17271634    }
    17281635
    17291636    // Allocate the bins, and initialize them to zero.  If there are N bounds,
    17301637    // then there are N-1 bins.
    1731     newHist->nums = psVectorAlloc( ( bounds->n ) - 1, PS_TYPE_U32 );
     1638    newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_U32);
    17321639    newHist->nums->n = newHist->nums->nalloc;
    1733     for ( i = 0;i < newHist->nums->n;i++ ) {
    1734         newHist->nums->data.U32[ i ] = 0;
     1640    for (i = 0; i < newHist->nums->n; i++) {
     1641        newHist->nums->data.U32[i] = 0;
    17351642    }
    17361643
    … …  
    17401647    newHist->uniform = false;
    17411648
    1742     return ( newHist );
    1743 }
    1744 
    1745 static void histogramFree( psHistogram *myHist )
    1746 {
    1747     psFree( myHist->bounds );
    1748     psFree( myHist->nums );
    1749 }
    1750 
     1649    return (newHist);
     1650}
     1651
     1652static void histogramFree(psHistogram * myHist)
     1653{
     1654    psFree(myHist->bounds);
     1655    psFree(myHist->nums);
     1656}
    17511657
    17521658/*****************************************************************************
    … …  
    17641670    The histogram structure "out".
    17651671 *****************************************************************************/
    1766 psHistogram *psVectorHistogram( psHistogram *out,
    1767                                 const psVector *restrict in,
    1768                                 const psVector *restrict mask,
    1769                                 unsigned int maskVal )
    1770 {
    1771     int i = 0;                                  // Loop index variable
    1772     int j = 0;                                  // Loop index variable
    1773     float binSize = 0.0;                        // Histogram bin size
    1774     int binNum = 0;                             // A temporary bin number
    1775     int numBins = 0;                            // The total number of bins
     1672psHistogram *psVectorHistogram(psHistogram * out,
     1673                               const psVector * restrict in,
     1674                               const psVector * restrict mask, unsigned int maskVal)
     1675{
     1676    int i = 0;                  // Loop index variable
     1677    int j = 0;                  // Loop index variable
     1678    float binSize = 0.0;        // Histogram bin size
     1679    int binNum = 0;             // A temporary bin number
     1680    int numBins = 0;            // The total number of bins
    17761681
    17771682    // NOTE: Verify that this is the correct action.
    1778     if ( out == NULL ) {
    1779         return ( NULL );
    1780     }
    1781 
     1683    if (out == NULL) {
     1684        return (NULL);
     1685    }
    17821686    // Check the specified output histogram for type psF32
    1783     if ( out->bounds->type.type != PS_TYPE_F32 ) {
    1784         psAbort( __func__,
    1785                  "Only data type PS_TYPE_F32 for the output.bounds member." );
    1786     }
    1787 
    1788     if ( out->nums->type.type != PS_TYPE_U32 ) {
    1789         psAbort( __func__,
    1790                  "Only data type PS_TYPE_U32 for output.nums member." );
    1791     }
    1792 
     1687    if (out->bounds->type.type != PS_TYPE_F32) {
     1688        psAbort(__func__, "Only data type PS_TYPE_F32 for the output.bounds member.");
     1689    }
     1690
     1691    if (out->nums->type.type != PS_TYPE_U32) {
     1692        psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member.");
     1693    }
    17931694    // NOTE: Verify that this is the correct action.
    1794     if ( in == NULL ) {
    1795         return ( out );
    1796     }
    1797 
    1798     if ( in->type.type != PS_TYPE_F32 ) {
    1799         psAbort( __func__,
    1800                  "Only data type PS_TYPE_F32 is currently supported (0x%x).",
    1801                  in->type.type );
    1802     }
    1803 
    1804     if ( mask != NULL ) {
    1805         if ( in->n != mask->n ) {
    1806             psAbort( __func__,
    1807                      "Vector data and vector mask are of different sizes." );
    1808         }
    1809         if ( mask->type.type != PS_TYPE_U8 ) {
    1810             psAbort( __func__, "Vector mask must be type PS_TYPE_U8" );
     1695    if (in == NULL) {
     1696        return (out);
     1697    }
     1698
     1699    if (in->type.type != PS_TYPE_F32) {
     1700        psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type);
     1701    }
     1702
     1703    if (mask != NULL) {
     1704        if (in->n != mask->n) {
     1705            psAbort(__func__, "Vector data and vector mask are of different sizes.");
     1706        }
     1707        if (mask->type.type != PS_TYPE_U8) {
     1708            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
    18111709        }
    18121710    }
    … …  
    18151713
    18161714    numBins = out->nums->n;
    1817     for ( i = 0;i < in->n;i++ ) {
     1715    for (i = 0; i < in->n; i++) {
    18181716        // Check if this pixel is masked, and if so, skip it.
    1819         if ( ( mask == NULL ) ||
    1820                 ( ( mask != NULL ) && ( !( mask->data.U8[ i ] & maskVal ) ) ) ) {
     1717        if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
    18211718            // Check if this pixel is below the minimum value, and if so
    18221719            // count it, then skip it.
    1823             if ( in->data.F32[ i ] < out->bounds->data.F32[ 0 ] ) {
     1720            if (in->data.F32[i] < out->bounds->data.F32[0]) {
    18241721                out->minNum++;
    18251722                // Check if this pixel is above the maximum value, and if so
    18261723                // count it, then skip it.
    1827             } else
    1828                 if ( in->data.F32[ i ] > out->bounds->data.F32[ numBins ] ) {
    1829                     out->maxNum++;
     1724            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
     1725                out->maxNum++;
     1726            } else {
     1727                // If this is a uniform histogram, determining the correct
     1728                // number is trivial.
     1729                if (out->uniform == true) {
     1730                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
     1731                    binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize);
     1732
     1733                    // NOTE: This next if-statement really shouldn't be necessary.
     1734                    // However, do to numerical lack of precision, we occasionally
     1735                    // produce a binNum outside the range of bins.
     1736                    if (binNum >= out->nums->n) {
     1737                        binNum = out->nums->n - 1;
     1738                    }
     1739
     1740                    (out->nums->data.S32[binNum])++;
     1741
     1742                    // If this is a non-uniform histogram, determining the correct
     1743                    // bin number requires a bit more work.
    18301744                } else {
    1831                     // If this is a uniform histogram, determining the correct
    1832                     // number is trivial.
    1833                     if ( out->uniform == true ) {
    1834                         binSize = out->bounds->data.F32[ 1 ] - out->bounds->data.F32[ 0 ];
    1835                         binNum = ( int ) ( ( in->data.F32[ i ] - out->bounds->data.F32[ 0 ] ) /
    1836                                            binSize );
    1837 
    1838                         // NOTE: This next if-statement really shouldn't be necessary.
    1839                         // However, do to numerical lack of precision, we occasionally
    1840                         // produce a binNum outside the range of bins.
    1841                         if ( binNum >= out->nums->n ) {
    1842                             binNum = out->nums->n - 1;
    1843                         }
    1844 
    1845                         ( out->nums->data.S32[ binNum ] ) ++;
    1846 
    1847                         // If this is a non-uniform histogram, determining the correct
    1848                         // bin number requires a bit more work.
    1849                     } else {
    1850                         // NOTE: This is slow.  Put a smarter algorithm here to
    1851                         // find the correct bin number (bin search, probably)
    1852                         for ( j = 0;j < ( out->bounds->n ) - 1;j++ ) {
    1853                             if ( ( out->bounds->data.S32[ j ] <= in->data.F32[ i ] ) &&
    1854                                     ( in->data.F32[ i ] <= out->bounds->data.S32[ j + 1 ] ) ) {
    1855                                 ( out->nums->data.S32[ j ] ) ++;
    1856                             }
     1745                    // NOTE: This is slow.  Put a smarter algorithm here to
     1746                    // find the correct bin number (bin search, probably)
     1747                    for (j = 0; j < (out->bounds->n) - 1; j++) {
     1748                        if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
     1749                                (in->data.F32[i] <= out->bounds->data.S32[j + 1])) {
     1750                            (out->nums->data.S32[j])++;
    18571751                        }
    18581752                    }
    18591753                }
    1860         }
    1861     }
    1862     return ( out );
     1754            }
     1755        }
     1756    }
     1757    return (out);
    18631758}
    18641759
    … …  
    18721767the various stat functions.
    18731768 *****************************************************************************/
    1874 psVector *p_psConvertToF32( psStats *stats,
    1875                             psVector *in,
    1876                             psVector *mask,
    1877                             unsigned int maskVal )
     1769psVector *p_psConvertToF32(psStats * stats, psVector * in, psVector * mask, unsigned int maskVal)
    18781770{
    18791771    int i = 0;
    18801772    psVector *tmp = NULL;
    18811773
    1882     if ( in->type.type == PS_TYPE_S32 ) {
    1883         tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
    1884         for ( i = 0;i < in->n;i++ ) {
    1885             tmp->data.F32[ i ] = ( float ) in->data.S32[ i ];
    1886         }
    1887     } else
    1888         if ( in->type.type == PS_TYPE_U32 ) {
    1889             tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
    1890             for ( i = 0;i < in->n;i++ ) {
    1891                 tmp->data.F32[ i ] = ( float ) in->data.U32[ i ];
    1892             }
    1893         } else
    1894             if ( in->type.type == PS_TYPE_F64 ) {
    1895                 tmp = psVectorAlloc( in->n, PS_TYPE_F32 );
    1896                 for ( i = 0;i < in->n;i++ ) {
    1897                     tmp->data.F32[ i ] = ( float ) in->data.F64[ i ];
    1898                 }
    1899             } else
    1900                 if ( in->type.type == PS_TYPE_F32 ) {
    1901                     // do nothing
    1902                 } else {
    1903                     psAbort( __func__, "unsupported vector type 0x%x\n", in->type.type );
    1904                 }
    1905     return ( tmp );
    1906 }
    1907 
     1774    if (in->type.type == PS_TYPE_S32) {
     1775        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     1776        for (i = 0; i < in->n; i++) {
     1777            tmp->data.F32[i] = (float)in->data.S32[i];
     1778        }
     1779    } else if (in->type.type == PS_TYPE_U32) {
     1780        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     1781        for (i = 0; i < in->n; i++) {
     1782            tmp->data.F32[i] = (float)in->data.U32[i];
     1783        }
     1784    } else if (in->type.type == PS_TYPE_F64) {
     1785        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     1786        for (i = 0; i < in->n; i++) {
     1787            tmp->data.F32[i] = (float)in->data.F64[i];
     1788        }
     1789    } else if (in->type.type == PS_TYPE_F32) {
     1790        // do nothing
     1791    } else {
     1792        psAbort(__func__, "unsupported vector type 0x%x\n", in->type.type);
     1793    }
     1794    return (tmp);
     1795}
    19081796
    19091797/******************************************************************************
    … …  
    19241812macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
    19251813 *****************************************************************************/
    1926 psStats *psVectorStats( psStats *stats,
    1927                         psVector *in,
    1928                         psVector *mask,
    1929                         unsigned int maskVal )
    1930 {
    1931     psVector * inF32;
     1814psStats *psVectorStats(psStats * stats, psVector * in, psVector * mask, unsigned int maskVal)
     1815{
     1816    psVector *inF32;
    19321817    int mustFreeTmp = 1;
    19331818
    19341819    // NOTE: Verify that this is the correct action.
    1935     if ( in == NULL ) {
    1936         return ( stats );
    1937     }
    1938     if ( stats == NULL ) {
    1939         return ( NULL );
    1940     }
    1941 
    1942     inF32 = p_psConvertToF32( stats, in, mask, maskVal );
    1943     if ( inF32 == NULL ) {
     1820    if (in == NULL) {
     1821        return (stats);
     1822    }
     1823    if (stats == NULL) {
     1824        return (NULL);
     1825    }
     1826
     1827    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
     1828    if (inF32 == NULL) {
    19441829        inF32 = in;
    19451830        mustFreeTmp = 0;
    19461831    }
    1947 
    19481832    // XXX: Should we abort if (stats->min == stats->max)?
    1949     if ( ( stats->options & PS_STAT_USE_RANGE ) &&
    1950             ( stats->min >= stats->max ) ) {
    1951         psAbort( __func__, "psVectorStats() called with range: %f to %f\n",
    1952                  stats->min, stats->max );
    1953     }
    1954 
    1955     //    PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
    1956     if ( mask != NULL ) {
    1957         PS_CHECK_NULL_VECTOR( mask );
    1958         PS_CHECK_EMPTY_VECTOR( mask );
    1959         PS_CHECK_VECTOR_SIZE_EQUAL( mask, in );
    1960         PS_CHECK_VECTOR_TYPE( mask, PS_TYPE_U8 );
    1961     }
    1962 
     1833    if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
     1834        psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
     1835    }
     1836    // PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
     1837    if (mask != NULL) {
     1838        PS_CHECK_NULL_VECTOR(mask);
     1839        PS_CHECK_EMPTY_VECTOR(mask);
     1840        PS_CHECK_VECTOR_SIZE_EQUAL(mask, in);
     1841        PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8);
     1842    }
    19631843    // ************************************************************************
    1964     if ( stats->options & PS_STAT_SAMPLE_MEAN ) {
    1965         p_psVectorSampleMean( in, mask, maskVal, stats );
    1966     }
    1967 
     1844    if (stats->options & PS_STAT_SAMPLE_MEAN) {
     1845        p_psVectorSampleMean(in, mask, maskVal, stats);
     1846    }
    19681847    // ************************************************************************
    1969     if ( stats->options & PS_STAT_SAMPLE_MEDIAN ) {
    1970         p_psVectorSampleMedian( in, mask, maskVal, stats );
    1971     }
    1972 
     1848    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
     1849        p_psVectorSampleMedian(in, mask, maskVal, stats);
     1850    }
    19731851    // ************************************************************************
    19741852    // NOTE: The Stdev calculation requires the mean.  Should we assume the
    1975     // mean has already been calculated?  Or should we always calculate it?
    1976     if ( stats->options & PS_STAT_SAMPLE_STDEV ) {
    1977         p_psVectorSampleMean( in, mask, maskVal, stats );
    1978         p_psVectorSampleStdev( in, mask, maskVal, stats );
    1979     }
    1980 
     1853    // mean has already been calculated? Or should we always calculate it?
     1854    if (stats->options & PS_STAT_SAMPLE_STDEV) {
     1855        p_psVectorSampleMean(in, mask, maskVal, stats);
     1856        p_psVectorSampleStdev(in, mask, maskVal, stats);
     1857    }
    19811858    // ************************************************************************
    1982     if ( stats->options & PS_STAT_SAMPLE_QUARTILE ) {
    1983         p_psVectorSampleQuartiles( in, mask, maskVal, stats );
    1984     }
    1985 
     1859    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
     1860        p_psVectorSampleQuartiles(in, mask, maskVal, stats);
     1861    }
    19861862    // Since the various robust stats quantities share much computation, they
    19871863    // are grouped together in a single private function:
    19881864    // p_psVectorRobustStats()
    1989     if ( ( stats->options & PS_STAT_ROBUST_MEAN ) ||
    1990             ( stats->options & PS_STAT_ROBUST_MEDIAN ) ||
    1991             ( stats->options & PS_STAT_ROBUST_MODE ) ||
    1992             ( stats->options & PS_STAT_ROBUST_STDEV ) ||
    1993             ( stats->options & PS_STAT_ROBUST_QUARTILE ) ) {
    1994         p_psVectorRobustStats( in, mask, maskVal, stats );
    1995     }
    1996 
    1997     if ( ( stats->options & PS_STAT_CLIPPED_MEAN ) ||
    1998             ( stats->options & PS_STAT_CLIPPED_STDEV ) ) {
    1999         p_psVectorClippedStats( in, mask, maskVal, stats );
    2000     }
    2001 
     1865    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
     1866            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
     1867            (stats->options & PS_STAT_ROBUST_MODE) ||
     1868            (stats->options & PS_STAT_ROBUST_STDEV) || (stats->options & PS_STAT_ROBUST_QUARTILE)) {
     1869        p_psVectorRobustStats(in, mask, maskVal, stats);
     1870    }
     1871
     1872    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
     1873        p_psVectorClippedStats(in, mask, maskVal, stats);
     1874    }
    20021875    // ************************************************************************
    2003     if ( stats->options & PS_STAT_MAX ) {
    2004         p_psVectorMax( in, mask, maskVal, stats );
    2005     }
    2006 
     1876    if (stats->options & PS_STAT_MAX) {
     1877        p_psVectorMax(in, mask, maskVal, stats);
     1878    }
    20071879    // ************************************************************************
    2008     if ( stats->options & PS_STAT_MIN ) {
    2009         p_psVectorMin( in, mask, maskVal, stats );
    2010     }
    2011 
    2012     if ( mustFreeTmp == 1 ) {
    2013         psFree( inF32 );
    2014     }
    2015     return ( stats );
    2016 }
     1880    if (stats->options & PS_STAT_MIN) {
     1881        p_psVectorMin(in, mask, maskVal, stats);
     1882    }
     1883
     1884    if (mustFreeTmp == 1) {
     1885        psFree(inF32);
     1886    }
     1887    return (stats);
     1888}
Note: See TracChangeset for help on using the changeset viewer.