IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1277


Ignore:
Timestamp:
Jul 22, 2004, 1:40:08 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psMinimize.c

    r1233 r1277  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-07-15 23:52:34 $
     11 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-07-22 23:40:08 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 */
    1616/*****************************************************************************/
    17 /* INCLUDE FILES            */
     17/* INCLUDE FILES                                                             */
    1818/*****************************************************************************/
    1919#include <stdlib.h>
    … …  
    4141#include "psMatrix.h"
    4242/*****************************************************************************/
    43 /* DEFINE STATEMENTS           */
     43/* DEFINE STATEMENTS                                                         */
    4444/*****************************************************************************/
    4545#define MAX_LMM_ITERATIONS 100
    … …  
    8585
    8686/*****************************************************************************/
    87 /* TYPE DEFINITIONS           */
     87/* TYPE DEFINITIONS                                                          */
    8888/*****************************************************************************/
    8989typedef struct
    … …  
    113113
    114114/*****************************************************************************/
    115 /* GLOBAL VARIABLES           */
     115/* GLOBAL VARIABLES                                                          */
    116116/*****************************************************************************/
    117117
    … …  
    119119
    120120/*****************************************************************************/
    121 /* FILE STATIC VARIABLES           */
     121/* FILE STATIC VARIABLES                                                     */
    122122/*****************************************************************************/
    123123
    … …  
    125125
    126126/*****************************************************************************/
    127 /* FUNCTION IMPLEMENTATION - LOCAL          */
     127/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    128128/*****************************************************************************/
    129129
    … …  
    558558
    559559/*****************************************************************************/
    560 /* FUNCTION IMPLEMENTATION - PUBLIC         */
     560/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    561561/*****************************************************************************/
    562562
    … …  
    567567 
    568568This routine must minimize an arbitrary function; it determines the set of
    569 parameters of that function such that the.
     569parameters of that function such that the ...
    570570 *****************************************************************************/
    571571psVector *
    … …  
    681681
    682682/******************************************************************************
    683     This routine must minimize an arbitrary function.
     683    This routine must determine the parameters of an arbitrary function
     684    such that they best fit the supplied data points.
    684685 *****************************************************************************/
    685686psVector *
    … …  
    795796    do {
    796797        iter++;
     798        for (i=0;i<initialGuess->n;i++) {
     799            printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
     800        }
    797801        // Perform an iteration of the GSL solver.
    798802        status = gsl_multifit_fdfsolver_iterate(s);
    799         //        printf ("psMinimize() status = %s\n", gsl_strerror(status));
     803        printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status));
     804        for (i=0;i<initialGuess->n;i++) {
     805            printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
     806        }
     807
    800808        // If there was a problem, abort.
    801809        if (status) {
    802             psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
     810            psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status));
    803811        }
    804812
    … …  
    811819        // as specified in the ADD.
    812820        *chiSq = gsl_blas_dnrm2(s->f);
     821        printf("psMinimize.c: chiSq is %.3f\n", *chiSq);
    813822        if (fabs(*chiSq - chiSqOld) < 1.0) {
    814823            status = GSL_SUCCESS;
    … …  
    855864    This routine must fit a polynomial of degree myPoly to the data points
    856865    (x, y) and return the coefficients of that polynomial, as well as the
    857     error for each data poiny (yErr).
     866    error for each data point (yErr).
    858867 
    859868NOTE: yErr is currently ignored.
    … …  
    877886    psVector *xSums = NULL;
    878887
     888    //    printf("psGetArrayPolynomial()\n");
     889    //    for (i=0;i<x->n;i++) {
     890    //        printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
     891    //    }
     892
    879893    PS_CHECK_NULL_1DPOLY(myPoly);
    880894    PS_CHECK_NULL_VECTOR(x);
    … …  
    929943    for(k=0;k<(polyOrder);k++) {
    930944        myPoly->coeff[k] = coeffs->data.F64[k];
    931     }
     945        //        printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
     946    }
     947
     948
     949    //    for (i=0;i<x->n;i++) {
     950    //        printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i], myPoly));
     951    //    }
    932952
    933953    psFree(A);
    … …  
    939959    psFree(xSums);
    940960
    941     return(NULL);
    942 }
     961    return(myPoly);
     962}

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

    r1233 r1277  
    55#include <float.h>
    66#include <math.h>
    7 
     7/*****************************************************************************/
     8/* INCLUDE FILES                                                             */
     9/*****************************************************************************/
    810#include "psMemory.h"
     11#include "psImage.h"
    912#include "psVector.h"
    1013#include "psTrace.h"
    … …  
    1215#include "psAbort.h"
    1316#include "psStats.h"
    14 
    15 #include "float.h"
    16 #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
    17 // will use robust statistical methods.
    18 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
    19 #define CLIPPED_NUM_ITER_LB 1
    20 #define CLIPPED_NUM_ITER_UB 10
    21 #define CLIPPED_SIGMA_LB 1.0
    22 #define CLIPPED_SIGMA_UB 10.0
    23 #define true 1
    24 #define false 0
    25 #define MYMAXFLOAT HUGE
    26 
    27 #ifndef DOXYGEN
     17#include "psSort.h"
     18#include "psMinimize.h"
     19#include "psFunctions.h"
     20=======
     21
     22    /*****************************************************************************/
     23    /* DEFINE STATEMENTS                                                         */
     24    /*****************************************************************************/
     25    #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
     26    // will use robust statistical methods.
     27    #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
     28    #define CLIPPED_NUM_ITER_LB 1
     29    #define CLIPPED_NUM_ITER_UB 10
     30    #define CLIPPED_SIGMA_LB 1.0
     31    #define CLIPPED_SIGMA_UB 10.0
     32    #define true 1
     33    #define false 0
     34    #define MYMAXFLOAT HUGE
     35    #define MAX_ITERATIONS 20
     36
     37    #ifndef DOXYGEN
     38    void p_psVectorRobustStats(const psVector *restrict myVector,
     39                               const psVector *restrict maskVector,
     40                               unsigned int maskVal,
     41                               psStats *stats);
     42#endif
     43
     44/** Preprocessor macro to generate error on an incorrect type */
     45#define PS_CHECK_VECTOR_TYPE(NAME, TYPE)                                                          \
     46if (NAME->type.type != TYPE) {                                                         \
     47    psError(__func__,"Invalid operation: %s has incorrect type.", #NAME);                          \
     48}
     49
     50/** Preprocessor macro to generate error on a NULL vector */
     51#define PS_CHECK_NULL_VECTOR(NAME)                                                          \
     52if (NAME == NULL || NAME->data.V == NULL) {                                                         \
     53    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     54}
     55
     56/** Preprocessor macro to generate error for zero length vector */
     57#define PS_CHECK_EMPTY_VECTOR(NAME)                                                          \
     58if (NAME->n < 1) {                                                                                  \
     59    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                             \
     60}
     61
     62/** Preprocessor macro to generate error on differing size vectors */
     63#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2)                                                          \
     64if (VEC1->n != VEC2->n) {               \
     65    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     66}
     67
     68/*****************************************************************************/
     69/* TYPE DEFINITIONS                                                          */
     70/*****************************************************************************/
     71
     72/*****************************************************************************/
     73/* GLOBAL VARIABLES                                                          */
     74/*****************************************************************************/
     75
     76// None
     77
     78/*****************************************************************************/
     79/* FILE STATIC VARIABLES                                                     */
     80/*****************************************************************************/
     81
     82// None
     83
     84/*****************************************************************************/
     85/* FUNCTION IMPLEMENTATION - LOCAL                                           */
     86/*****************************************************************************/
     87
     88/*****************************************************************************
     89p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
     90function that simply prints a vector to STDOUT.  Used primarily for
     91debugging.
     92 *****************************************************************************/
     93
     94void p_psVectorPrint(psVector *myVector,
     95                     psVector *maskVector,
     96                     unsigned int maskVal,
     97                     psStats *stats)
     98{
     99    int i = 0;                                  // Loop index variable.
     100
     101    for (i=0;i<myVector->n;i++) {
     102        if (maskVector != NULL)
     103            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
     104        else
     105            printf("Element %d is %f\n", i, myVector->data.F32[i]);
     106    }
     107}
     108
     109/******************************************************************************
     110 ******************************************************************************
     111 ******************************************************************************
     112    MISC PRIVATE STATISTICAL FUNCTIONS
     113 
     114    NOTE: it is assumed that any call to these statistical functions will
     115    have been preceded by a call to the psVectorStats() function.  Various
     116    sanity tests will only be performed in psVectorStats().
     117        Is the mask vector the same length as the data vector?
     118        Is the mask vector of type PS_TYPE_U8?
     119        Is the stats data structure NULL?
     120        Is the in data structure NULL?
     121        Is the in data structure of type PS_TYPE_F32?
     122 ******************************************************************************
     123 ******************************************************************************
     124 *****************************************************************************/
     125
     126
     127/******************************************************************************
     128p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
     129mean of the input vector.
     130Inputs
     131Returns
     132    NULL
     133ASSUMPTION: the mean is always calculated exactly.  Robust means are never
     134calculated in this routine.
     135 *****************************************************************************/
     136void p_psVectorSampleMean(const psVector *restrict myVector,
     137                          const psVector *restrict maskVector,
     138                          unsigned int maskVal,
     139                          psStats *stats)
     140{
     141    int i = 0;                                  // Loop index variable
     142    float mean = 0.0;                           // The mean
     143    int count = 0;                              // # of points in this mean?
     144    float rangeMin = 0.0;                       // Exclude data below this
     145    float rangeMax = 0.0;                       // Exclude date above this
     146
     147    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
     148    // loop.
     149    if (stats->options & PS_STAT_USE_RANGE) {
     150        rangeMin = stats->min;
     151        rangeMax = stats->max;
     152        if (maskVector != NULL) {
     153            for (i=0;i<myVector->n;i++) {
     154                // Check if the data is with the specified range
     155                if (!(maskVal & maskVector->data.U8[i]) &&
     156                        (rangeMin <= myVector->data.F32[i]) &&
     157                        (myVector->data.F32[i] <= rangeMax)) {
     158                    mean+= myVector->data.F32[i];
     159                    count++;
     160                }
     161            }
     162            mean/= (float) count;
     163        } else {
     164            for (i=0;i<myVector->n;i++) {
     165                if ((rangeMin <= myVector->data.F32[i]) &&
     166                        (myVector->data.F32[i] <= rangeMax)) {
     167                    mean+= myVector->data.F32[i];
     168                    count++;
     169                }
     170            }
     171            mean/= (float) count;
     172        }
     173    } else {
     174        if (maskVector != NULL) {
     175            for (i=0;i<myVector->n;i++) {
     176                if (!(maskVal & maskVector->data.U8[i])) {
     177                    mean+= myVector->data.F32[i];
     178                    count++;
     179                }
     180            }
     181            mean/= (float) count;
     182        } else {
     183            for (i=0;i<myVector->n;i++) {
     184                mean+= myVector->data.F32[i];
     185            }
     186            mean/= (float) myVector->n;
     187        }
     188    }
     189
     190    stats->sampleMean = mean;
     191}
     192
     193/******************************************************************************
     194p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     195max of the input vector.
     196Inputs
     197    myVector
     198    maskVector
     199    maskVal
     200    stats
     201Returns
     202    NULL
     203 *****************************************************************************/
     204void p_psVectorMax(const psVector *restrict myVector,
     205                   const psVector *restrict maskVector,
     206                   unsigned int maskVal,
     207                   psStats *stats)
     208{
     209    int i = 0;                                  // Loop index variable
     210    float max = -MYMAXFLOAT;                    // The calculated maximum
     211    float rangeMin = 0.0;                       // Exclude data below this
     212    float rangeMax = 0.0;                       // Exclude date above this
     213
     214    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
     215    if (stats->options & PS_STAT_USE_RANGE) {
     216        rangeMin = stats->min;
     217        rangeMax = stats->max;
     218        if (maskVector != NULL) {
     219            for (i=0;i<myVector->n;i++) {
     220                if (!(maskVal & maskVector->data.U8[i])) {
     221                    if ((myVector->data.F32[i] > max) &&
     222                            (rangeMin <= myVector->data.F32[i]) &&
     223                            (myVector->data.F32[i] <= rangeMax)) {
     224                        max = myVector->data.F32[i];
     225                    }
     226                }
     227            }
     228        } else {
     229            for (i=0;i<myVector->n;i++) {
     230                if ((myVector->data.F32[i] > max) &&
     231                        (rangeMin <= myVector->data.F32[i]) &&
     232                        (myVector->data.F32[i] <= rangeMax)) {
     233                    max = myVector->data.F32[i];
     234                }
     235            }
     236        }
     237    } else {
     238        if (maskVector != NULL) {
     239            for (i=0;i<myVector->n;i++) {
     240                if (!(maskVal & maskVector->data.U8[i])) {
     241                    if (myVector->data.F32[i] > max) {
     242                        max = myVector->data.F32[i];
     243                    }
     244                }
     245            }
     246        } else {
     247            for (i=0;i<myVector->n;i++) {
     248                if (myVector->data.F32[i] > max) {
     249                    max = myVector->data.F32[i];
     250                }
     251            }
     252        }
     253    }
     254
     255    stats->max = max;
     256}
     257
     258/******************************************************************************
     259p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
     260minimum of the input vector.
     261Inputs
     262    myVector
     263    maskVector
     264    maskVal
     265    stats
     266Returns
     267    NULL
     268 *****************************************************************************/
     269void p_psVectorMin(const psVector *restrict myVector,
     270                   const psVector *restrict maskVector,
     271                   unsigned int maskVal,
     272                   psStats *stats)
     273{
     274    int i = 0;                                  // Loop index variable
     275    float min = MYMAXFLOAT;                     // The calculated maximum
     276    float rangeMin = 0.0;                       // Exclude data below this
     277    float rangeMax = 0.0;                       // Exclude date above this
     278
     279    if (stats->options & PS_STAT_USE_RANGE) {
     280        rangeMin = stats->min;
     281        rangeMax = stats->max;
     282        if (maskVector != NULL) {
     283            for (i=0;i<myVector->n;i++) {
     284                if (!(maskVal & maskVector->data.U8[i])) {
     285                    if ((myVector->data.F32[i] < min) &&
     286                            (rangeMin <= myVector->data.F32[i]) &&
     287                            (myVector->data.F32[i] <= rangeMax)) {
     288                        min = myVector->data.F32[i];
     289                    }
     290                }
     291            }
     292        } else {
     293            for (i=0;i<myVector->n;i++) {
     294                if ((myVector->data.F32[i] < min) &&
     295                        (rangeMin <= myVector->data.F32[i]) &&
     296                        (myVector->data.F32[i] <= rangeMax)) {
     297                    min = myVector->data.F32[i];
     298                }
     299            }
     300        }
     301    } else {
     302        if (maskVector != NULL) {
     303            for (i=0;i<myVector->n;i++) {
     304                if (!(maskVal & maskVector->data.U8[i])) {
     305                    if (myVector->data.F32[i] < min) {
     306                        min = myVector->data.F32[i];
     307                    }
     308                }
     309            }
     310        } else {
     311            for (i=0;i<myVector->n;i++) {
     312                if (myVector->data.F32[i] < min) {
     313                    min = myVector->data.F32[i];
     314                }
     315            }
     316        }
     317    }
     318
     319    stats->min = min;
     320}
     321
     322/******************************************************************************
     323p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
     324number of non-masked pixels in the vector that fall within the min/max
     325range, if given.
     326Inputs
     327    myVector
     328    maskVector
     329    maskVal
     330    stats
     331Returns
     332    NULL
     333 *****************************************************************************/
     334int p_psVectorNValues(const psVector *restrict myVector,
     335                      const psVector *restrict maskVector,
     336                      unsigned int maskVal,
     337                      psStats *stats)
     338{
     339    int i = 0;                                  // Loop index variable
     340    int numData = 0;                            // The number of data points
     341    float rangeMin = 0.0;                       // Exclude data below this
     342    float rangeMax = 0.0;                       // Exclude date above this
     343
     344    if (stats->options & PS_STAT_USE_RANGE) {
     345        rangeMin = stats->min;
     346        rangeMax = stats->max;
     347        if (maskVector != NULL) {
     348            for (i=0;i<myVector->n;i++) {
     349                if (!(maskVal & maskVector->data.U8[i]) &&
     350                        (rangeMin <= myVector->data.F32[i]) &&
     351                        (myVector->data.F32[i] <= rangeMax)) {
     352                    numData++;
     353                }
     354            }
     355        } else {
     356            for (i=0;i<myVector->n;i++) {
     357                if ((rangeMin <= myVector->data.F32[i]) &&
     358                        (myVector->data.F32[i] <= rangeMax)) {
     359                    numData++;
     360                }
     361            }
     362        }
     363    } else {
     364        rangeMin = stats->min;
     365        rangeMax = stats->max;
     366        if (maskVector != NULL) {
     367            for (i=0;i<myVector->n;i++) {
     368                if (!(maskVal & maskVector->data.U8[i])) {
     369                    numData++;
     370                }
     371            }
     372        } else {
     373            numData = myVector->n;
     374        }
     375    }
     376    return(numData);
     377}
     378
     379
     380
     381/******************************************************************************
     382p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
     383median of the input vector.
     384Inputs
     385    myVector
     386    maskVector
     387    maskVal
     388    stats
     389Returns
     390    NULL
     391 *****************************************************************************/
     392void p_psVectorSampleMedian(const psVector *restrict myVector,
     393                            const psVector *restrict maskVector,
     394                            unsigned int maskVal,
     395                            psStats *stats)
     396{
     397    psVector *unsortedVector = NULL;            // Temporary vector
     398    psVector *sortedVector = NULL;              // Temporary vector
     399    int i = 0;                                  // Loop index variable
     400    int count = 0;                              // # of points in this mean?
     401    int nValues = 0;                            // # of points in vector
     402    float rangeMin = 0.0;                       // Exclude data below this
     403    float rangeMax = 0.0;                       // Exclude date above this
     404    psStats *stats2 = NULL;                     // Temporary stats structure
     405
     406
     407    // Determine if the number of data points exceed a threshold which will
     408    // cause to generate robust stats, as opposed to exact stats.
     409
     410    if (myVector->n > stats->sampleLimit) {
     411        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     412
     413        // Calculate the robust quartiles.
     414        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     415        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     416
     417        // Store the robust quartiles into the sample quartile members.
     418        stats->sampleMedian = stats2->robustMedian;
     419
     420        // Free temporary data buffers.
     421        psFree(stats2);
     422
     423        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     424        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     425
     426        return;
     427    }
     428
     429    // Determine how many data points fit inside this min/max range
     430    // and are not masked, IF the maskVector is not NULL>
     431    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
     432
     433    // Allocate temporary vectors for the data.
     434    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     435    unsortedVector->n = unsortedVector->nalloc;
     436
     437    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     438    sortedVector->n = sortedVector->nalloc;
     439
     440    // Determine if we must only use data points within a min/max range.
     441    if (stats->options & PS_STAT_USE_RANGE) {
     442        rangeMin = stats->min;
     443        rangeMax = stats->max;
     444
     445        // Store all non-masked data points within the min/max range
     446        // into the temporary vectors.
     447        count = 0;
     448        if (maskVector != NULL) {
     449            for (i=0;i<myVector->n;i++) {
     450                if (!(maskVal & maskVector->data.U8[i]) &&
     451                        (rangeMin <= myVector->data.F32[i]) &&
     452                        (myVector->data.F32[i] <= rangeMax)) {
     453                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
     454                }
     455            }
     456        } else {
     457            for (i=0;i<myVector->n;i++) {
     458                if ((rangeMin <= myVector->data.F32[i]) &&
     459                        (myVector->data.F32[i] <= rangeMax)) {
     460                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     461                }
     462            }
     463        }
     464    } else {
     465        // Store all non-masked data points into the temporary vectors.
     466        count = 0;
     467        if (maskVector != NULL) {
     468            for (i=0;i<myVector->n;i++) {
     469                if (!(maskVal & maskVector->data.U8[i])) {
     470                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     471                }
     472            }
     473        } else {
     474            for (i=0;i<myVector->n;i++) {
     475                unsortedVector->data.F32[i] = myVector->data.F32[i];
     476            }
     477        }
     478    }
     479    // Sort the temporary vectors.
     480    psVectorSort(sortedVector, unsortedVector);
     481
     482    // Calculate the median exactly.
     483    if (0 == (nValues % 2)) {
     484        stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
     485                                     sortedVector->data.F32[nValues/2]);
     486    } else {
     487        stats->sampleMedian = sortedVector->data.F32[nValues/2];
     488    }
     489
     490    // Free the temporary data structures.
     491    psFree(unsortedVector);
     492    psFree(sortedVector);
     493}
     494
     495/******************************************************************************
     496    This routine smoothes the data in the input robustHistogram with a
     497    Gaussian of width sigma.
     498 *****************************************************************************/
     499psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
     500                                       float sigma)
     501{
     502    int i = 0;                                  // Loop index variable
     503    int j = 0;                                  // Loop index variable
     504    float denom = 0.0;                           // Temporary variable
     505    float expo = 0.0;                            // Temporary variable
     506    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
     507    psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
     508
     509    for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
     510        if (fabs(sigma) >= FLT_EPSILON) {
     511            // If sigma does not equal zero, then we use Gaussian smoothing.
     512            #ifdef  DARWIN
     513            denom = (float) sqrt(2.0 * M_PI * sigma * sigma);
     514            #else
     515
     516            denom = sqrtf(2.0 * M_PI * sigma * sigma);
     517            #endif
     518
     519            expo = - (float) ((i-GAUSS_WIDTH) * (i-GAUSS_WIDTH));
     520            expo/= (2.0 * sigma * sigma);
     521            gaussianCoefs[i] = exp(expo/denom);
     522
     523            // NOTE: Gaussian smoothing just isn't working with low sigma
     524            // values.  The problem is that the Gaussian coefficients are
     525            // all zero, except for the middle coefficient, which is exactly
     526            // one.  Therefore, I'm using boxcar smoothing.
     527            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float) GAUSS_WIDTH));
     528            //            printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
     529        } else {
     530            /* If sigma equals zero (all pixels have the same value)
     531             * the above code will divide by zero.  Therefore, we don't need
     532             * to smooth the data.
     533             */
     534            for(i=0;i<robustHistogram->nums->n;i++) {
     535                smooth->data.F32[i] = (float) robustHistogram->nums->data.S32[i];
     536            }
     537            return(smooth);
     538        }
     539    }
     540
     541    for(i=0;i<robustHistogram->nums->n;i++) {
     542        smooth->data.F32[i] = 0.0;
     543        for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
     544            if (((j+i) >= 0) && ((j+i) < smooth->n)) {
     545                smooth->data.F32[i]+= (gaussianCoefs[j+GAUSS_WIDTH] *
     546                                       (float) robustHistogram->nums->data.S32[j+i]);
     547            }
     548        }
     549    }
     550    return(smooth);
     551}
     552
     553/******************************************************************************
     554p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
     555the upper and/or lower quartiles of the input vector.
     556Inputs
     557    myVector
     558    maskVector
     559    maskVal
     560    stats
     561Returns
     562    NULL
     563 *****************************************************************************/
     564void p_psVectorSampleQuartiles(const psVector *restrict myVector,
     565                               const psVector *restrict maskVector,
     566                               unsigned int maskVal,
     567                               psStats *stats)
     568{
     569    psVector *unsortedVector = NULL;            // Temporary vector
     570    psVector *sortedVector = NULL;              // Temporary vector
     571    int i = 0;                                  // Loop index variable
     572    int count = 0;                              // # of points in this mean?
     573    int nValues = 0;                            // # data points
     574    float rangeMin = 0.0;                       // Exclude data below this
     575    float rangeMax = 0.0;                       // Exclude date above this
     576
     577    // Determine if the number of data points exceed a threshold which will
     578    // cause to generate robust stats, as opposed to exact stats.
     579    if (myVector->n > stats->sampleLimit) {
     580        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     581        psStats *stats2 = NULL;
     582        // Calculate the robust quartiles.
     583        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
     584        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     585
     586        // Store the robust quartiles into the sample quartile members.
     587        stats->sampleUQ = stats2->robustUQ;
     588        stats->sampleLQ = stats2->robustLQ;
     589
     590        // Free temporary data buffers.
     591        psFree(stats2);
     592
     593        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     594        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     595
     596        return;
     597    }
     598
     599    // Determine how many data points fit inside this min/max range
     600    // and are not maxed, IF the maskVector is not NULL>
     601    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
     602
     603    // Allocate temporary vectors for the data.
     604    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     605    unsortedVector->n = unsortedVector->nalloc;
     606    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     607    sortedVector->n = sortedVector->nalloc;
     608
     609    // Determine if we must only use data points within a min/max range.
     610    if (stats->options & PS_STAT_USE_RANGE) {
     611        rangeMin = stats->min;
     612        rangeMax = stats->max;
     613        // Store all non-masked data points within the min/max range
     614        // into the temporary vectors.
     615        count = 0;
     616        if (maskVector != NULL) {
     617            for (i=0;i<myVector->n;i++) {
     618                if (!(maskVal & maskVector->data.U8[i]) &&
     619                        (rangeMin <= myVector->data.F32[i]) &&
     620                        (myVector->data.F32[i] <= rangeMax)) {
     621                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     622                }
     623            }
     624        } else {
     625            for (i=0;i<myVector->n;i++) {
     626                if ((rangeMin <= myVector->data.F32[i]) &&
     627                        (myVector->data.F32[i] <= rangeMax)) {
     628                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     629                }
     630            }
     631        }
     632    } else {
     633        // Store all non-masked data points into the temporary vectors.
     634        count = 0;
     635        if (maskVector != NULL) {
     636            for (i=0;i<myVector->n;i++) {
     637                if (!(maskVal & maskVector->data.U8[i])) {
     638                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     639                }
     640            }
     641        } else {
     642            for (i=0;i<myVector->n;i++) {
     643                unsortedVector->data.F32[i] = myVector->data.F32[i];
     644            }
     645        }
     646    }
     647
     648    // Sort the temporary vectors.
     649    psVectorSort(sortedVector, unsortedVector);
     650
     651    // Calculate the quartile points exactly.
     652    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     653    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
     654
     655    // Free the temporary data structures.
     656    psFree(unsortedVector);
     657    psFree(sortedVector);
     658    // NOTE: This is the
     659}
     660
     661
     662/******************************************************************************
     663p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
     664stdev of the input vector.
     665Inputs
     666    myVector
     667    maskVector
     668    maskVal
     669    stats
     670Returns
     671    NULL
     672 
     673NOTE: the mean is always calculated exactly.  Robust means are never
     674calculated in this routine.
     675 *****************************************************************************/
     676void p_psVectorSampleStdev(const psVector *restrict myVector,
     677                           const psVector *restrict maskVector,
     678                           unsigned int maskVal,
     679                           psStats *stats)
     680{
     681    int i = 0;                                  // Loop index variable
     682    int countInt = 0;                           // # of data points being used
     683    float countFloat = 0.0;                     // # of data points being used
     684    float mean = 0.0;                           // The mean
     685    float diff = 0.0;                           // Used in calculating stdev
     686    float sumSquares = 0.0;                     // temporary variable
     687    float sumDiffs = 0.0;                       // temporary variable
     688    float rangeMin = 0.0;                       // Exclude data below this
     689    float rangeMax = 0.0;                       // Exclude date above this
     690
     691    // This procedure requires the mean.  If it has not been already
     692    // calculated, then call p_psVectorSampleMean()
     693    if (0 != isnan(stats->sampleMean)) {
     694        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     695    }
     696    mean = stats->sampleMean;
     697
     698    if (stats->options & PS_STAT_USE_RANGE) {
     699        if (maskVector != NULL) {
     700            for (i=0;i<myVector->n;i++) {
     701                if (!(maskVal & maskVector->data.U8[i]) &&
     702                        (rangeMin <= myVector->data.F32[i]) &&
     703                        (myVector->data.F32[i] <= rangeMax)) {
     704                    diff = myVector->data.F32[i] - mean;
     705                    sumSquares+= (diff * diff);
     706                    sumDiffs+= diff;
     707                    countInt++;
     708                }
     709            }
     710        } else {
     711            for (i=0;i<myVector->n;i++) {
     712                if ((rangeMin <= myVector->data.F32[i]) &&
     713                        (myVector->data.F32[i] <= rangeMax)) {
     714                    diff = myVector->data.F32[i] - mean;
     715                    sumSquares+= (diff * diff);
     716                    sumDiffs+= diff;
     717                    countInt++;
     718                }
     719            }
     720            countInt = myVector->n;
     721        }
     722    } else {
     723        if (maskVector != NULL) {
     724            for (i=0;i<myVector->n;i++) {
     725                if (!(maskVal & maskVector->data.U8[i])) {
     726                    diff = myVector->data.F32[i] - mean;
     727                    sumSquares+= (diff * diff);
     728                    sumDiffs+= diff;
     729                    countInt++;
     730                }
     731            }
     732        } else {
     733            for (i=0;i<myVector->n;i++) {
     734                diff = myVector->data.F32[i] - mean;
     735                sumSquares+= (diff * diff);
     736                sumDiffs+= diff;
     737                countInt++;
     738            }
     739            countInt = myVector->n;
     740        }
     741    }
     742    countFloat = (float) countInt;
     743
     744    #ifdef DARWIN
     745
     746    stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
     747                                        sumDiffs/countFloat))/ (countFloat-1));
     748    #else
     749
     750    stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
     751                                 sumDiffs/countFloat))/ (countFloat-1));
     752    #endif
     753}
     754
     755/******************************************************************************
     756p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
     757clipped stats (mean or stdev) of the input vector.
     758 
     759Inputs
     760    myVector
     761    maskVector
     762    maskVal
     763    stats
     764Returns
     765    NULL
     766 *****************************************************************************/
     767void p_psVectorClippedStats(const psVector *restrict myVector,
     768                            const psVector *restrict maskVector,
     769                            unsigned int maskVal,
     770                            psStats *stats)
     771{
     772    int i = 0;                                  // Loop index variable
     773    int j = 0;                                  // Loop index variable
     774    float clippedMean = 0.0;                    // self-explanatory
     775    float clippedStdev = 0.0;                   // self-explanatory
     776    float oldStanMean = 0.0;                    // Temporary variable
     777    float oldStanStdev = 0.0;                   // Temporary variable
     778    psVector *tmpMask = NULL;                   // Temporary vector
     779
     780    // Endure that stats->clipIter is within the proper range.
     781    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
     782            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
     783        psAbort(__func__, "Unallowed value for clipIter (%d).\n",
     784                stats->clipIter);
     785    }
     786
     787    // Endure that stats->clipSigma is within the proper range.
     788    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
     789            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
     790        psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
     791                stats->clipSigma);
     792    }
     793
     794    // We allocate a temporary mask vector since during the iterative
     795    // steps that follow, we will be masking off additional data points.
     796    // However, we do no want to modify the original mask vector.
     797    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
     798    tmpMask->n = myVector->n;
     799
     800    // If we were called with a mask vector, then initialize the temporary
     801    // mask vector with those values.
     802    if (maskVector != NULL) {
     803        for (i=0;i<tmpMask->n;i++) {
     804            tmpMask->data.U8[i] = maskVector->data.U8[i];
     805        }
     806    }
     807
     808    // 1. Compute the sample median.
     809    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
     810    // median here, not the mean.
     811    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
     812
     813    // 2. Compute the sample standard deviation.
     814    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
     815
     816    // 3. Use the sample median as the first estimator of the mean X.
     817    clippedMean = stats->sampleMean;
     818
     819    // 4. Use the sample stdev as the first estimator of the mean stdev.
     820    clippedStdev = stats->sampleStdev;
     821
     822    // Must save the old sampleMean and sampleStdev since the following code
     823    // block overwrites them.
     824    oldStanMean = stats->sampleMean;
     825    oldStanStdev = stats->sampleStdev;
     826
     827    // 5. Repeat N times:
     828    for (i=0;i<stats->clipIter;i++) {
     829        for (j=0;j<myVector->n;j++) {
     830            // a) Exclude all values x_i for which |x_i - x| > K * stdev
     831            if (fabs(myVector->data.F32[j] - clippedMean) >
     832                    (stats->clipSigma * clippedStdev)) {
     833                tmpMask->data.U8[i] = 0xff;
     834            }
     835            // b) compute new mean and stdev
     836            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
     837            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
     838
     839            // c) Use the new mean for x
     840            clippedMean = stats->sampleMean;
     841
     842            // d) Use the new stdev for stdev
     843            clippedStdev = stats->sampleStdev;
     844        }
     845
     846    }
     847    stats->sampleMean = oldStanMean;
     848    stats->sampleStdev= oldStanStdev;
     849
     850    // 7. The last calcuated value of x is the cliped mean.
     851    if (stats->options & PS_STAT_CLIPPED_MEAN) {
     852        stats->clippedMean = clippedMean;
     853    }
     854
     855    // 8. The last calcuated value of stdev is the cliped stdev.
     856    if (stats->options & PS_STAT_CLIPPED_STDEV) {
     857        stats->clippedStdev = clippedStdev;
     858    }
     859
     860    psFree(tmpMask);
     861}
     862
     863/*****************************************************************************
     864 *****************************************************************************/
     865void p_psNormalizeVector(psVector *myData)
     866{
     867    float min = (float) HUGE;
     868    float max = (float) -HUGE;
     869    float range = 0.0;
     870    int i = 0;
     871
     872    for (i=0;i<myData->n;i++) {
     873        if (myData->data.F32[i] < min) {
     874            min = myData->data.F32[i];
     875        }
     876        if (myData->data.F32[i] > max) {
     877            max = myData->data.F32[i];
     878        }
     879    }
     880
     881    range = max - min;
     882    for (i=0;i<myData->n;i++) {
     883        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
     884    }
     885}
     886
     887
     888float p_psGaussian(const psVector *restrict myData,
     889                   const psVector *restrict myParams)
     890{
     891    float x = myData->data.F32[0];
     892    float mean = myParams->data.F32[0];
     893    float stdev = myParams->data.F32[1];
     894    float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev));
     895    tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
     896
     897    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     898    return(tmp);
     899}
     900
     901float p_psGaussianDeriv(const psVector *restrict myData,
     902                        const psVector *restrict myParams,
     903                        int whichParam)
     904{
     905    float x = myData->data.F32[0];
     906    float mean = myParams->data.F32[0];
     907    float stdev = myParams->data.F32[1];
     908    float tmp = 0.0;
     909
     910    if (whichParam == 0) {
     911        // Return the derivative w.r.t. the mean.
     912        tmp = (x - mean) * p_psGaussian(myData, myParams);
     913        tmp/= (stdev * stdev);
     914    } else if (whichParam == 1) {
     915        // Return the derivative w.r.t. the stdev.
     916        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
     917        tmp/= (stdev * stdev * stdev);
     918    }
     919    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
     920
     921    return(tmp);
     922}
     923
     924
     925float p_psQuadratic(const psVector *restrict myParams,
     926                    const psVector *restrict myCoords)
     927{
     928    float x = myCoords->data.F32[0];
     929    float A = myParams->data.F32[0];
     930    float B = myParams->data.F32[1];
     931    float C = myParams->data.F32[2];
     932    float tmp = 0.0;
     933
     934    tmp = (A * x * x) + (B * x) + C;
     935    return(tmp);
     936}
     937
     938float p_psQuadraticDeriv(const psVector *restrict myParams,
     939                         const psVector *restrict myCoords,
     940                         int whichParamDeriv)
     941{
     942    float x = myCoords->data.F32[0];
     943    float tmp = 0.0;
     944
     945    if (whichParamDeriv == 0) {
     946        tmp = x * x;
     947    } else if (whichParamDeriv == 1) {
     948        tmp = x;
     949    } else if (whichParamDeriv == 2) {
     950        tmp = 1.0;
     951    }
     952
     953    return(tmp);
     954}
     955
     956/******************************************************************************
     957p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
     958as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
     959polynomials here) and a range of x-values for which it is defined:
     960[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
     961that f(x) == midpoint.  This functions uses a binary-search algorithm on the
     962range and assumes that the polnomial is monotonically increasing within that
     963range.
     964 *****************************************************************************/
     965float p_ps1DPolyMedian(psPolynomial1D *myPoly,
     966                       float rangeLow,
     967                       float rangeHigh,
     968                       float getThisValue)
     969{
     970    int numIterations=0;
     971    float midpoint = 0.0;
     972    float oldMidpoint = 1.0;
     973    float f = 0.0;
     974
     975    //  printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
     976
     977    while (numIterations < MAX_ITERATIONS) {
     978        midpoint = (rangeHigh + rangeLow) / 2.0;
     979        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
     980            return(midpoint);
     981        }
     982        oldMidpoint = midpoint;
     983
     984        f = psEvalPolynomial1D(midpoint, myPoly);
     985        //        printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
     986        if (fabs(f - getThisValue) <= FLT_EPSILON) {
     987            return(midpoint);
     988        }
     989
     990        if (f > getThisValue) {
     991            rangeHigh = midpoint;
     992        } else {
     993            rangeLow = midpoint;
     994        }
     995        numIterations++;
     996    }
     997    return(midpoint);
     998}
     999
     1000/******************************************************************************
     1001p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
     1002as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
     1003polynomials here) and a range of x-values for which it is defined:
     1004[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
     1005that f(x) == midpoint.  This functions uses a binary-search algorithm on the
     1006range and assumes that the polnomial is monotonically increasing within that
     1007range.
     1008 *****************************************************************************/
     1009float p_psFitQuadratic(psHistogram *robustHistogram,
     1010                       int binNum,
     1011                       float fitFloat)
     1012{
     1013    /*
     1014        if ((binNum > 0) &&
     1015            (binNum < (robustHistogram->nums->n+1))) {
     1016            x->data.F64[0] = (double) 0.5 *
     1017                      (robustHistogram->bounds->data.F32[binNum-1] +
     1018                       robustHistogram->bounds->data.F32[binNum]);
     1019            x->data.F64[1] = (double) 0.5 *
     1020                      (robustHistogram->bounds->data.F32[binNum] +
     1021                       robustHistogram->bounds->data.F32[binNum+1]);
     1022            x->data.F64[2] = (double) 0.5 *
     1023                      (robustHistogram->bounds->data.F32[binNum+1] +
     1024                       robustHistogram->bounds->data.F32[binNum+2]);
     1025     
     1026            y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1];
     1027            y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum];
     1028            y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1];
     1029     
     1030            if (!((y->data.F64[0] <= fitFloat) &&
     1031                 (fitFloat <= y->data.F64[2]))) {
     1032                psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
     1033            }
     1034     
     1035            yErr->data.F64[0] = 1.0;
     1036            yErr->data.F64[1] = 1.0;
     1037            yErr->data.F64[2] = 1.0;
     1038     
     1039            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1040            return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat);
     1041         } else {
     1042            return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] +
     1043                          robustHistogram->bounds->data.F32[binNum]));
     1044        }
     1045    */
     1046    return(0.0);
     1047}
     1048
     1049#define PS_PRINT_VECTOR(NAME) \
     1050for (int my_i=0;my_i<NAME->n;my_i++) { \
     1051    printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \
     1052} \
     1053printf("\n"); \
     1054
     1055/******************************************************************************
     1056p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
     1057calculates a variety of robust stat measures:
     1058PS_STAT_ROBUST_MEAN
     1059PS_STAT_ROBUST_MEDIAN
     1060PS_STAT_ROBUST_MODE
     1061PS_STAT_ROBUST_STDEV
     1062PS_STAT_ROBUST_QUARTILE
     1063I have included all that computation in a single function, as opposed to
     1064breaking it across several functions for one primary reason:
     1065they all
     1066require the same basic initial processing steps (calculate the histogram,
     1067etc.)
     1068however there is no currently defined means, in the SDRS, to
     1069specify this computation as a separate step.  If the above robust stat
     1070measures were calcualted in separate functiosn, then much of the initial
     1071processing would be duplicated.
     1072Inputs
     1073myVector
     1074maskVector
     1075maskVal
     1076stats
     1077Returns
     1078NULL
     1079*****************************************************************************/
    281080void p_psVectorRobustStats(const psVector *restrict myVector,
    291081                           const psVector *restrict maskVector,
    301082                           unsigned int maskVal,
    31                            psStats *stats);
    32 #endif
     1083                           psStats *stats)
     1084{
     1085    psHistogram *robustHistogram = NULL;
     1086    psVector *robustHistogramVector = NULL;
     1087    float binSize = 0.0;                        // Size of the histogram bins
     1088    int LQBinNum = -1;                          // Bin num for lower quartile
     1089    int UQBinNum = -1;                          // Bin num for upper quartile
     1090    int i = 0;                                  // Loop index variable
     1091    int maxBinNum = 0;
     1092    float maxBinCount = 0.0;
     1093    float dL = 0.0;
     1094    int numBins = 0;
     1095    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
     1096    //    psImage *domain;
     1097    //    psVector *errors;
     1098    //    psVector *data;
     1099    //    psVector *initialGuess;
     1100    //    psVector *theParams;
     1101    //    float chiSq=0.0;
     1102    //    float max = -HUGE;
     1103    //    float max_pos;
     1104    float myMean = 0.0;
     1105    float myStdev = 0.0;
     1106    float countFloat = 0.0;
     1107    float diff = 0.0;
     1108    float sumSquares = 0.0;
     1109    float sumDiffs = 0.0;
     1110    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1111    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1112    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1113    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
     1114    psVector *cumulativeRobustSumsFullRange = NULL;
     1115    psVector *cumulativeRobustSumsDlRange = NULL;
     1116    float sumRobust = 0.0;
     1117    float sumN50 = 0.0;
     1118    float sumNfit = 0.0;
     1119    float cumulativeMedian = 0.0;
     1120
     1121    // Compute the initial bin size of the robust histogram.  This is done
     1122    // by computing the clipped standard deviation of the vector, and dividing
     1123    // that by 10.0;
     1124    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
     1125    binSize = tmpStats->clippedStdev / 10.0f;
     1126
     1127    // If stats->clippedStdev == 0.0, then all data elements have the same
     1128    // value.  Therefore, we can set the appropiate results and return.
     1129    if (fabs(binSize) <= FLT_EPSILON) {
     1130        if (stats->options & PS_STAT_ROBUST_MEAN) {
     1131            stats->robustMean = stats->clippedMean;
     1132        }
     1133        if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
     1134            stats->robustMedian = stats->clippedMean;
     1135        }
     1136        if  (stats->options & PS_STAT_ROBUST_MODE) {
     1137            stats->robustMode = stats->clippedMean;
     1138        }
     1139        if  (stats->options & PS_STAT_ROBUST_STDEV) {
     1140            stats->robustStdev = 0.0;
     1141        }
     1142        if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
     1143            stats->robustUQ = stats->clippedMean;
     1144            stats->robustLQ = stats->clippedMean;
     1145        }
     1146        // XXX: Set these to the number of unmasked data points?
     1147        stats->robustNfit = 0.0;
     1148        stats->robustN50 = 0.0;
     1149        psFree(tmpStats);
     1150        return;
     1151    }
     1152
     1153    // Determine minimum and maximum values in the data vector.
     1154    if (isnan(stats->min)) {
     1155        p_psVectorMin(myVector, maskVector, maskVal, stats);
     1156    }
     1157    if (isnan(stats->max)) {
     1158        p_psVectorMax(myVector, maskVector, maskVal, stats);
     1159    }
     1160
     1161    // Create the histogram structure.  NOTE: we can not specify the bin size
     1162    // precisely since the argument to psHistogramAlloc() is the number of
     1163    // bins, not the binSize.  Also, if we get here, we know that
     1164    // binSize != 0.0.
     1165    numBins = (int) ((stats->max - stats->min) / binSize);
     1166    robustHistogram = psHistogramAlloc(stats->min,
     1167                                       stats->max,
     1168                                       numBins);
     1169
     1170    // Populate the histogram array.
     1171    psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
     1172
     1173    // Smooth the histogram.
     1174    // XXX: is that the right stdev?
     1175    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
     1176                            tmpStats->clippedStdev/4.0f);
     1177
     1178    // The following was necessary to fit a gaussian to the data, since
     1179    // gaussian functions produce data between 0.0 and 1.0.
     1180    // p_psNormalizeVector(robustHistogramVector);
     1181
     1182    /**************************************************************************
     1183    Determine the lower/upper quartiles.
     1184    **************************************************************************/
     1185    // We define a vector called "cumulativeRobustSums..." where the value at
     1186    // index position i is equal to the sum of bins 0:i.  This will be used
     1187    // now and later in determining the lower/upper quartiles.
     1188    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1189    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
     1190    for (i=1;i<robustHistogramVector->n;i++) {
     1191        cumulativeRobustSumsFullRange->data.F32[i] =
     1192            cumulativeRobustSumsFullRange->data.F32[i-1] +
     1193            robustHistogramVector->data.F32[i];
     1194    }
     1195    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1];
     1196
     1197    // Determine the bin number containing the lower quartile point.
     1198    LQBinNum = -1;
     1199    for (i=0;i<cumulativeRobustSumsFullRange->n;i++) {
     1200        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust/4.0)) {
     1201            LQBinNum = i;
     1202            break;
     1203        }
     1204    }
     1205
     1206    // Determine the bin number containing the upper quartile point.
     1207    UQBinNum = -1;
     1208    for (i=cumulativeRobustSumsFullRange->n-1;i>=0;i--) {
     1209        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0*sumRobust/4.0)) {
     1210            UQBinNum = i;
     1211            break;
     1212        }
     1213    }
     1214
     1215    if ((LQBinNum == -1) ||
     1216            (UQBinNum == -1)) {
     1217        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
     1218    }
     1219    /**************************************************************************
     1220    Determine the mode in the range LQ:UQ.
     1221    **************************************************************************/
     1222    // Determine the bin with the peak value in the range LQ to UQ.
     1223    maxBinNum = LQBinNum;
     1224    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
     1225    sumN50 = (float) robustHistogram->nums->data.S32[LQBinNum];
     1226    for (i=LQBinNum+1;i<=UQBinNum;i++) {
     1227        if (robustHistogramVector->data.F32[i] > maxBinCount) {
     1228            maxBinNum = i;
     1229            maxBinCount = robustHistogramVector->data.F32[i];
     1230        }
     1231        sumN50+= (float) robustHistogram->nums->data.S32[i];
     1232    }
     1233    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
     1234    dL = (UQBinNum - LQBinNum) / 4;
     1235
     1236    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
     1237
     1238    /**************************************************************************
     1239    Determine the mean/stdev for the bins in the range mode-dL to mode+dL
     1240    **************************************************************************/
     1241    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1242    for (i=0;i<robustHistogramVector->n;i++) {
     1243        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
     1244    }
     1245    sumNfit = 0.0;
     1246    cumulativeMedian = 0.0;
     1247    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1248        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1249            cumulativeRobustSumsDlRange->data.F32[i] =
     1250                cumulativeRobustSumsDlRange->data.F32[i-1] +
     1251                robustHistogramVector->data.F32[i];
     1252            cumulativeMedian+= robustHistogramVector->data.F32[i];
     1253            sumNfit+= (float) robustHistogramVector->data.S32[i];
     1254        }
     1255    }
     1256
     1257    // Calculate the mean of the smoothed robust histogram in the range
     1258    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     1259    // that bin (this is a non-exact approximation).
     1260    myMean = 0.0;
     1261    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1262        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1263            myMean+= (robustHistogramVector->data.F32[i]) * 0.5 *
     1264                     (robustHistogram->bounds->data.F32[i+1] +
     1265                      robustHistogram->bounds->data.F32[i]);
     1266            countFloat+= robustHistogramVector->data.F32[i];
     1267        }
     1268    }
     1269    myMean/= countFloat;
     1270
     1271    // Calculate the stdev of the smoothed robust histogram in the range
     1272    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     1273    // that bin.
     1274    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1275        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1276            diff = (0.5 * (robustHistogram->bounds->data.F32[i+1] +
     1277                           robustHistogram->bounds->data.F32[i])) - myMean;
     1278            sumSquares+= diff * diff * robustHistogramVector->data.F32[i];
     1279            sumDiffs+= diff * robustHistogramVector->data.F32[i];
     1280        }
     1281    }
     1282    myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1));
     1283
     1284    /**************************************************************************
     1285    Set the appropriate members in the output stats struct.
     1286    **************************************************************************/
     1287    if (stats->options & PS_STAT_ROBUST_MEAN) {
     1288        stats->robustMean = myMean;
     1289    }
     1290
     1291    if  (stats->options & PS_STAT_ROBUST_MODE) {
     1292        stats->robustMode = 0.5 *
     1293                            (robustHistogram->bounds->data.F32[maxBinNum] +
     1294                             robustHistogram->bounds->data.F32[maxBinNum+1]);
     1295    }
     1296
     1297    if  (stats->options & PS_STAT_ROBUST_STDEV) {
     1298        stats->robustStdev = myStdev;
     1299    }
     1300
     1301    // To determine the median (and later, the lower/upper quartile), we fit
     1302    // a quadratic to the three bins surrounding the bin containing the median.
     1303    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
     1304    // the cumulative number of data points in all bins up to, and including,
     1305    // this bin.  We then solve the quadratic for
     1306
     1307    if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
     1308        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n-1))) {
     1309            x->data.F64[0] = (double) 0.5 *
     1310                             (robustHistogram->bounds->data.F32[maxBinNum-1] +
     1311                              robustHistogram->bounds->data.F32[maxBinNum]);
     1312            x->data.F64[1] = (double) 0.5 *
     1313                             (robustHistogram->bounds->data.F32[maxBinNum] +
     1314                              robustHistogram->bounds->data.F32[maxBinNum+1]);
     1315            x->data.F64[2] = (double) 0.5 *
     1316                             (robustHistogram->bounds->data.F32[maxBinNum+1] +
     1317                              robustHistogram->bounds->data.F32[maxBinNum+2]);
     1318
     1319            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum-1];
     1320            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
     1321            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum+1];
     1322
     1323            // Ensure that cumulativeMedian/2 is actually within the range of the bins
     1324            // we are using.
     1325            cumulativeMedian*= 0.5;
     1326            if (!((y->data.F64[0] <= cumulativeMedian) &&
     1327                    (cumulativeMedian <= y->data.F64[2]))) {
     1328                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
     1329                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
     1330            }
     1331            // XXX: yErr is not currently used by psGetArrayPolynomial().  We
     1332            // may have to set this meaningfully later.
     1333            yErr->data.F64[0] = 1.0;
     1334            yErr->data.F64[1] = 1.0;
     1335            yErr->data.F64[2] = 1.0;
     1336
     1337            // Determine the coefficients of the polynomial.
     1338            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1339            // Call p_ps1DPolyMedian(), which does a binary search on the
     1340            // polynomial, looking for the value x such that
     1341            // f(x) = cumulativeMedian.
     1342            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0],
     1343                                                   x->data.F64[2], cumulativeMedian);
     1344        } else {
     1345            // If the mode is the first/last histogram bin, then simply use
     1346            // the midpoint of that bin.
     1347            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum+1] +
     1348                                         robustHistogram->bounds->data.F32[maxBinNum]);
     1349        }
     1350    }
     1351
     1352    // The lower/upper quartile calculations are very similar to the median
     1353    // calculations.  We fit a quadratic to the array containing the
     1354    // cumulative data points in each bin, then search for x such that
     1355    // f(x) equals the lower/upper quartile exactly.
     1356    //
     1357    if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
     1358        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1];
     1359
     1360        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n-1))) {
     1361            x->data.F64[0] = (double) 0.5 *
     1362                             (robustHistogram->bounds->data.F32[LQBinNum-1] +
     1363                              robustHistogram->bounds->data.F32[LQBinNum]);
     1364            x->data.F64[1] = (double) 0.5 *
     1365                             (robustHistogram->bounds->data.F32[LQBinNum] +
     1366                              robustHistogram->bounds->data.F32[LQBinNum+1]);
     1367            x->data.F64[2] = (double) 0.5 *
     1368                             (robustHistogram->bounds->data.F32[LQBinNum+1] +
     1369                              robustHistogram->bounds->data.F32[LQBinNum+2]);
     1370
     1371            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum-1];
     1372            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
     1373            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum+1];
     1374
     1375            if (!((y->data.F64[0] <= (countFloat/4.0)) &&
     1376                    ((countFloat/4.0) <= y->data.F64[2]))) {
     1377                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
     1378            }
     1379
     1380            yErr->data.F64[0] = 1.0;
     1381            yErr->data.F64[1] = 1.0;
     1382            yErr->data.F64[2] = 1.0;
     1383
     1384            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1385            stats->robustLQ = p_ps1DPolyMedian(myPoly,
     1386                                               x->data.F64[0],
     1387                                               x->data.F64[2],
     1388                                               countFloat/4.0);
     1389
     1390        } else {
     1391            // If the LQ is the first/last histogram bin, then simply use
     1392            // the midpoint of that bin.
     1393            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum+1] +
     1394                                     robustHistogram->bounds->data.F32[LQBinNum]);
     1395        }
     1396
     1397        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n-1))) {
     1398            x->data.F64[0] = (double) 0.5 *
     1399                             (robustHistogram->bounds->data.F32[UQBinNum-1] +
     1400                              robustHistogram->bounds->data.F32[UQBinNum]);
     1401            x->data.F64[1] = (double) 0.5 *
     1402                             (robustHistogram->bounds->data.F32[UQBinNum] +
     1403                              robustHistogram->bounds->data.F32[UQBinNum+1]);
     1404            x->data.F64[2] = (double) 0.5 *
     1405                             (robustHistogram->bounds->data.F32[UQBinNum+1] +
     1406                              robustHistogram->bounds->data.F32[UQBinNum+2]);
     1407
     1408            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum-1];
     1409            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
     1410            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum+1];
     1411
     1412            if (!((y->data.F64[0] <= (3.0 * countFloat/4.0)) &&
     1413                    ((3.0 * countFloat/4.0) <= y->data.F64[2]))) {
     1414                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
     1415            }
     1416
     1417            yErr->data.F64[0] = 1.0;
     1418            yErr->data.F64[1] = 1.0;
     1419            yErr->data.F64[2] = 1.0;
     1420
     1421            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1422            stats->robustUQ = p_ps1DPolyMedian(myPoly,
     1423                                               x->data.F64[0],
     1424                                               x->data.F64[2],
     1425                                               3.0*countFloat/4.0);
     1426        } else {
     1427            // If the UQ is the first/last histogram bin, then simply use
     1428            // the midpoint of that bin.
     1429            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum+1] +
     1430                                     robustHistogram->bounds->data.F32[UQBinNum]);
     1431        }
     1432    }
     1433    stats->robustNfit = sumNfit;
     1434    stats->robustN50 = sumN50;
     1435
     1436    psFree(x);
     1437    psFree(y);
     1438    psFree(yErr);
     1439    psFree(tmpStats);
     1440    psFree(robustHistogram);
     1441    psFree(myPoly);
     1442    psFree(cumulativeRobustSumsFullRange);
     1443    psFree(cumulativeRobustSumsDlRange);
     1444}
     1445
     1446
     1447
     1448/*
     1449void p_ps_FitTheGaussian()
     1450{
     1451 
     1452    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     1453    // NOTE: This code uses the psMinimize.c functions to perform
     1454    // the fit.  It doesn't quite work, 100% of the time.  For the time being
     1455    // I am commenting this code out, and replacing it with code which
     1456    // calculates the mean directly on the robustHistogram.
     1457 
     1458    domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
     1459    errors = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1460    data = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1461    initialGuess = psVectorAlloc(2, PS_TYPE_F32);
     1462 
     1463    max = -HUGE;
     1464    max_pos = -1;
     1465    for (i=0;i<robustHistogramVector->n;i++) {
     1466        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
     1467                                        robustHistogram->bounds->data.F32[i]);
     1468        data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
     1469        errors->data.F32[i] = 1.0;
     1470        //printf("DATA (%.2f, %.2f)\n", domain->data.F32[i][0], data->data.F32[i]);
     1471        if (data->data.F32[i] > max) {
     1472            max = data->data.F32[i];
     1473            max_pos = domain->data.F32[i][0];
     1474        }
     1475    }
     1476 
     1477    initialGuess->data.F32[0] = max_pos;
     1478    initialGuess->data.F32[1] = 1.0;
     1479    printf("Initial (mean. stdev) is (%.3f, %.3f)\n", initialGuess->data.F32[0], initialGuess->data.F32[1]);
     1480 
     1481    printf("Calling psMinimizeChi2()\n");
     1482    theParams = psMinimizeChi2(p_psGaussian,
     1483                               p_psGaussianDeriv,
     1484                               domain,
     1485                               data,
     1486                               errors,
     1487                               initialGuess,
     1488                               NULL,
     1489                               &chiSq);
     1490    printf("Called psMinimizeChi2()\n");
     1491    printf("p_psVectorRobustStats() is (%f, %f)\n", theParams->data.F32[0], theParams->data.F32[1]);
     1492}
     1493*/
     1494
     1495
     1496/*****************************************************************************/
     1497/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     1498/*****************************************************************************/
    331499
    341500static void histogramFree(psHistogram *myHist);
    … …  
    921558        return(NULL);
    931559    }
     1560    if (n < 0) {
     1561        psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n);
     1562    }
    941563
    951564    // NOTE: Verify that this is the correct action.
    … …  
    1071576    // Calculate the bounds for each bin.
    1081577    binSize = (upper - lower) / (float) n;
     1578    // NOTE: Is the following necessary?  It prevents the max data point
     1579    // from being in a non-existant bin.
     1580    binSize+= FLT_EPSILON;
    1091581    for (i=0;i<n+1;i++) {
    1101582        newHist->bounds->data.F32[i] = lower + (binSize * (float) i);
    … …  
    1271599
    1281600/******************************************************************************
    129 psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure
    130 with the specifed bounds.  The number of elements in the bounds vector is n.
    131 Therefore, the number of bins is n-1.
     1601psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure
     1602with the specifed bounds.
     1603 
    1321604Input:
    1331605    bounds
    … …  
    2381710
    2391711    numBins = out->nums->n;
    240 
    2411712    for (i=0;i<in->n;i++) {
    2421713        // Check if this pixel is masked, and if so, skip it.
    … …  
    2561727                if (out->uniform == true) {
    2571728                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    258 
    2591729                    binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    2601730                                    binSize);
     1731
     1732                    // NOTE: This next if-statement really shouldn't be necessary.
     1733                    // However, do to numerical lack of precision, we occasionally
     1734                    // produce a binNum outside the range of bins.
     1735                    if (binNum >= out->nums->n) {
     1736                        binNum = out->nums->n-1;
     1737                    }
     1738
    2611739                    (out->nums->data.S32[binNum])++;
     1740
    2621741                    // If this is a non-uniform histogram, determining the correct
    2631742                    // bin number requires a bit more work.
    … …  
    2771756    return(out);
    2781757}
    279 
    280 /*****************************************************************************
    281 p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
    282 function that simply prints a vector to STDOUT.  Used primarily for
    283 debugging.
    284  *****************************************************************************/
    285 
    286 void p_psVectorPrint(psVector *myVector,
    287                      psVector *maskVector,
    288                      unsigned int maskVal,
    289                      psStats *stats)
    290 {
    291     int i = 0;                                  // Loop index variable.
    292 
    293     for (i=0;i<myVector->n;i++) {
    294         if (maskVector != NULL)
    295             printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
    296         else
    297             printf("Element %d is %f\n", i, myVector->data.F32[i]);
    298     }
    299 }
    300 
    301 
    302 /******************************************************************************
    303  ******************************************************************************
    304  ******************************************************************************
    305     MISC PRIVATE STATISTICAL FUNCTIONS
    306  
    307     NOTE: it is assumed that any call to these statistical functions will
    308     have been preceded by a call to the psVectorStats() function.  Various
    309     sanity tests will only be performed in psVectorStats().
    310         Is the mask vector the same length as the data vector?
    311         Is the mask vector of type PS_TYPE_U8?
    312         Is the stats data structure NULL?
    313         Is the in data structure NULL?
    314         Is the in data structure of type PS_TYPE_F32?
    315  ******************************************************************************
    316  ******************************************************************************
    317  *****************************************************************************/
    318 
    319 
    320 /******************************************************************************
    321 p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
    322 mean of the input vector.
    323 Inputs
    324 Returns
    325     NULL
    326 ASSUMPTION: the mean is always calculated exactly.  Robust means are never
    327 calculated in this routine.
    328  *****************************************************************************/
    329 void p_psVectorSampleMean(const psVector *restrict myVector,
    330                           const psVector *restrict maskVector,
    331                           unsigned int maskVal,
    332                           psStats *stats)
    333 {
    334     int i = 0;                                  // Loop index variable
    335     float mean = 0.0;                           // The mean
    336     int count = 0;                              // # of points in this mean?
    337     float rangeMin = 0.0;                       // Exclude data below this
    338     float rangeMax = 0.0;                       // Exclude date above this
    339 
    340     // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
    341     // loop.
    342     if (stats->options & PS_STAT_USE_RANGE) {
    343         rangeMin = stats->min;
    344         rangeMax = stats->max;
    345         if (maskVector != NULL) {
    346             for (i=0;i<myVector->n;i++) {
    347                 // Check if the data is with the specified range
    348                 if (!(maskVal & maskVector->data.U8[i]) &&
    349                         (rangeMin <= myVector->data.F32[i]) &&
    350                         (myVector->data.F32[i] <= rangeMax)) {
    351                     mean+= myVector->data.F32[i];
    352                     count++;
    353                 }
    354             }
    355             mean/= (float) count;
    356         } else {
    357             for (i=0;i<myVector->n;i++) {
    358                 if ((rangeMin <= myVector->data.F32[i]) &&
    359                         (myVector->data.F32[i] <= rangeMax)) {
    360                     mean+= myVector->data.F32[i];
    361                     count++;
    362                 }
    363             }
    364             mean/= (float) count;
    365         }
    366     } else {
    367         if (maskVector != NULL) {
    368             for (i=0;i<myVector->n;i++) {
    369                 if (!(maskVal & maskVector->data.U8[i])) {
    370                     mean+= myVector->data.F32[i];
    371                     count++;
    372                 }
    373             }
    374             mean/= (float) count;
    375         } else {
    376             for (i=0;i<myVector->n;i++) {
    377                 mean+= myVector->data.F32[i];
    378             }
    379             mean/= (float) myVector->n;
    380         }
    381     }
    382 
    383     stats->sampleMean = mean;
    384 }
    385 
    386 /******************************************************************************
    387 p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
    388 max of the input vector.
    389 Inputs
    390     myVector
    391     maskVector
    392     maskVal
    393     stats
    394 Returns
    395     NULL
    396  *****************************************************************************/
    397 void p_psVectorMax(const psVector *restrict myVector,
    398                    const psVector *restrict maskVector,
    399                    unsigned int maskVal,
    400                    psStats *stats)
    401 {
    402     int i = 0;                                  // Loop index variable
    403     float max = -MYMAXFLOAT;                    // The calculated maximum
    404     float rangeMin = 0.0;                       // Exclude data below this
    405     float rangeMax = 0.0;                       // Exclude date above this
    406 
    407     // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    408     if (stats->options & PS_STAT_USE_RANGE) {
    409         rangeMin = stats->min;
    410         rangeMax = stats->max;
    411         if (maskVector != NULL) {
    412             for (i=0;i<myVector->n;i++) {
    413                 if (!(maskVal & maskVector->data.U8[i])) {
    414                     if ((myVector->data.F32[i] > max) &&
    415                             (rangeMin <= myVector->data.F32[i]) &&
    416                             (myVector->data.F32[i] <= rangeMax)) {
    417                         max = myVector->data.F32[i];
    418                     }
    419                 }
    420             }
    421         } else {
    422             for (i=0;i<myVector->n;i++) {
    423                 if ((myVector->data.F32[i] > max) &&
    424                         (rangeMin <= myVector->data.F32[i]) &&
    425                         (myVector->data.F32[i] <= rangeMax)) {
    426                     max = myVector->data.F32[i];
    427                 }
    428             }
    429         }
    430     } else {
    431         if (maskVector != NULL) {
    432             for (i=0;i<myVector->n;i++) {
    433                 if (!(maskVal & maskVector->data.U8[i])) {
    434                     if (myVector->data.F32[i] > max) {
    435                         max = myVector->data.F32[i];
    436                     }
    437                 }
    438             }
    439         } else {
    440             for (i=0;i<myVector->n;i++) {
    441                 if (myVector->data.F32[i] > max) {
    442                     max = myVector->data.F32[i];
    443                 }
    444             }
    445         }
    446     }
    447 
    448     stats->max = max;
    449 }
    450 
    451 /******************************************************************************
    452 p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
    453 minimum of the input vector.
    454 Inputs
    455     myVector
    456     maskVector
    457     maskVal
    458     stats
    459 Returns
    460     NULL
    461  *****************************************************************************/
    462 void p_psVectorMin(const psVector *restrict myVector,
    463                    const psVector *restrict maskVector,
    464                    unsigned int maskVal,
    465                    psStats *stats)
    466 {
    467     int i = 0;                                  // Loop index variable
    468     float min = MYMAXFLOAT;                     // The calculated maximum
    469     float rangeMin = 0.0;                       // Exclude data below this
    470     float rangeMax = 0.0;                       // Exclude date above this
    471 
    472     if (stats->options & PS_STAT_USE_RANGE) {
    473         rangeMin = stats->min;
    474         rangeMax = stats->max;
    475         if (maskVector != NULL) {
    476             for (i=0;i<myVector->n;i++) {
    477                 if (!(maskVal & maskVector->data.U8[i])) {
    478                     if ((myVector->data.F32[i] < min) &&
    479                             (rangeMin <= myVector->data.F32[i]) &&
    480                             (myVector->data.F32[i] <= rangeMax)) {
    481                         min = myVector->data.F32[i];
    482                     }
    483                 }
    484             }
    485         } else {
    486             for (i=0;i<myVector->n;i++) {
    487                 if ((myVector->data.F32[i] < min) &&
    488                         (rangeMin <= myVector->data.F32[i]) &&
    489                         (myVector->data.F32[i] <= rangeMax)) {
    490                     min = myVector->data.F32[i];
    491                 }
    492             }
    493         }
    494     } else {
    495         if (maskVector != NULL) {
    496             for (i=0;i<myVector->n;i++) {
    497                 if (!(maskVal & maskVector->data.U8[i])) {
    498                     if (myVector->data.F32[i] < min) {
    499                         min = myVector->data.F32[i];
    500                     }
    501                 }
    502             }
    503         } else {
    504             for (i=0;i<myVector->n;i++) {
    505                 if (myVector->data.F32[i] < min) {
    506                     min = myVector->data.F32[i];
    507                 }
    508             }
    509         }
    510     }
    511 
    512     stats->min = min;
    513 }
    514 
    515 /******************************************************************************
    516 p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
    517 number of non-masked pixels in the vector that fall within the min/max
    518 range, if given.
    519 Inputs
    520     myVector
    521     maskVector
    522     maskVal
    523     stats
    524 Returns
    525     NULL
    526  *****************************************************************************/
    527 int p_psVectorNValues(const psVector *restrict myVector,
    528                       const psVector *restrict maskVector,
    529                       unsigned int maskVal,
    530                       psStats *stats)
    531 {
    532     int i = 0;                                  // Loop index variable
    533     int numData = 0;                            // The number of data points
    534     float rangeMin = 0.0;                       // Exclude data below this
    535     float rangeMax = 0.0;                       // Exclude date above this
    536 
    537     if (stats->options & PS_STAT_USE_RANGE) {
    538         rangeMin = stats->min;
    539         rangeMax = stats->max;
    540         if (maskVector != NULL) {
    541             for (i=0;i<myVector->n;i++) {
    542                 if (!(maskVal & maskVector->data.U8[i]) &&
    543                         (rangeMin <= myVector->data.F32[i]) &&
    544                         (myVector->data.F32[i] <= rangeMax)) {
    545                     numData++;
    546                 }
    547             }
    548         } else {
    549             for (i=0;i<myVector->n;i++) {
    550                 if ((rangeMin <= myVector->data.F32[i]) &&
    551                         (myVector->data.F32[i] <= rangeMax)) {
    552                     numData++;
    553                 }
    554             }
    555         }
    556     } else {
    557         rangeMin = stats->min;
    558         rangeMax = stats->max;
    559         if (maskVector != NULL) {
    560             for (i=0;i<myVector->n;i++) {
    561                 if (!(maskVal & maskVector->data.U8[i])) {
    562                     numData++;
    563                 }
    564             }
    565         } else {
    566             numData = myVector->n;
    567         }
    568     }
    569     return(numData);
    570 }
    571 
    572 
    573 
    574 /******************************************************************************
    575 p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
    576 median of the input vector.
    577 Inputs
    578     myVector
    579     maskVector
    580     maskVal
    581     stats
    582 Returns
    583     NULL
    584  *****************************************************************************/
    585 void p_psVectorSampleMedian(const psVector *restrict myVector,
    586                             const psVector *restrict maskVector,
    587                             unsigned int maskVal,
    588                             psStats *stats)
    589 {
    590     psVector *unsortedVector = NULL;            // Temporary vector
    591     psVector *sortedVector = NULL;              // Temporary vector
    592     int i = 0;                                  // Loop index variable
    593     int count = 0;                              // # of points in this mean?
    594     int nValues = 0;                            // # of points in vector
    595     float rangeMin = 0.0;                       // Exclude data below this
    596     float rangeMax = 0.0;                       // Exclude date above this
    597     psStats *stats2 = NULL;                     // Temporary stats structure
    598 
    599 
    600     // Determine if the number of data points exceed a threshold which will
    601     // cause to generate robust stats, as opposed to exact stats.
    602 
    603     if (myVector->n > stats->sampleLimit) {
    604         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    605 
    606         // Calculate the robust quartiles.
    607         stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    608         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    609 
    610         // Store the robust quartiles into the sample quartile members.
    611         stats->sampleMedian = stats2->robustMedian;
    612 
    613         // Free temporary data buffers.
    614         psFree(stats2);
    615 
    616         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    617         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    618 
    619         return;
    620     }
    621 
    622     // Determine how many data points fit inside this min/max range
    623     // and are not masked, IF the maskVector is not NULL>
    624     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    625 
    626     // Allocate temporary vectors for the data.
    627     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    628     unsortedVector->n = unsortedVector->nalloc;
    629 
    630     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
    631     sortedVector->n = sortedVector->nalloc;
    632 
    633     // Determine if we must only use data points within a min/max range.
    634     if (stats->options & PS_STAT_USE_RANGE) {
    635         rangeMin = stats->min;
    636         rangeMax = stats->max;
    637 
    638         // Store all non-masked data points within the min/max range
    639         // into the temporary vectors.
    640         count = 0;
    641         if (maskVector != NULL) {
    642             for (i=0;i<myVector->n;i++) {
    643                 if (!(maskVal & maskVector->data.U8[i]) &&
    644                         (rangeMin <= myVector->data.F32[i]) &&
    645                         (myVector->data.F32[i] <= rangeMax)) {
    646                     unsortedVector->data.F32[count++] = maskVector->data.F32[i];
    647                 }
    648             }
    649         } else {
    650             for (i=0;i<myVector->n;i++) {
    651                 if ((rangeMin <= myVector->data.F32[i]) &&
    652                         (myVector->data.F32[i] <= rangeMax)) {
    653                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    654                 }
    655             }
    656         }
    657     } else {
    658         // Store all non-masked data points into the temporary vectors.
    659         count = 0;
    660         if (maskVector != NULL) {
    661             for (i=0;i<myVector->n;i++) {
    662                 if (!(maskVal & maskVector->data.U8[i])) {
    663                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    664                 }
    665             }
    666         } else {
    667             for (i=0;i<myVector->n;i++) {
    668                 unsortedVector->data.F32[i] = myVector->data.F32[i];
    669             }
    670         }
    671     }
    672     // Sort the temporary vectors.
    673     psVectorSort(sortedVector, unsortedVector);
    674 
    675     // Calculate the median exactly.
    676     if (0 == (nValues % 2)) {
    677         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
    678                                      sortedVector->data.F32[nValues/2]);
    679     } else {
    680         stats->sampleMedian = sortedVector->data.F32[nValues/2];
    681     }
    682 
    683     // Free the temporary data structures.
    684     psFree(unsortedVector);
    685     psFree(sortedVector);
    686 }
    687 
    688 /******************************************************************************
    689     This routine smoothes the data in the input robustHistogram with a
    690     Gaussian of width sigma.
    691  *****************************************************************************/
    692 void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
    693                                   float sigma)
    694 {
    695     int i = 0;                                  // Loop index variable
    696     int j = 0;                                  // Loop index variable
    697     float tmpf = 0.0;                           // Temporary variable
    698     float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
    699 
    700     for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
    701         if (fabs(sigma) <= FLT_EPSILON) {
    702             // If sigma does not equal zero, then we use Gaussian smoothing.
    703             #ifdef  DARWIN
    704             tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma);
    705             #else
    706 
    707             tmpf = sqrtf(2.0f * M_PI * sigma * sigma);
    708             #endif
    709 
    710             gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) *
    711                                              ((float) (i-GAUSS_WIDTH))) /
    712                                             (2.0f * sigma * sigma)) / tmpf;
    713         } else {
    714             /* If sigma equals zero (all pixels have the same value)
    715              * the above code will divide by zero.  Therefore, we don't need
    716              * to smooth the data.
    717              */
    718             return;
    719         }
    720     }
    721 
    722     for(i=0;i<robustHistogram->nums->n;i++) {
    723         for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
    724             if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) {
    725                 robustHistogram->nums->data.S32[j+i]+=
    726                     (gaussianCoefs[j+GAUSS_WIDTH] *
    727                      (float) robustHistogram->nums->data.S32[j+i]);
    728             }
    729         }
    730     }
    731 }
    732 
    733 /******************************************************************************
    734 p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
    735 the upper and/or lower quartiles of the input vector.
    736 Inputs
    737     myVector
    738     maskVector
    739     maskVal
    740     stats
    741 Returns
    742     NULL
    743  *****************************************************************************/
    744 void p_psVectorSampleQuartiles(const psVector *restrict myVector,
    745                                const psVector *restrict maskVector,
    746                                unsigned int maskVal,
    747                                psStats *stats)
    748 {
    749     psVector *unsortedVector = NULL;            // Temporary vector
    750     psVector *sortedVector = NULL;              // Temporary vector
    751     int i = 0;                                  // Loop index variable
    752     int count = 0;                              // # of points in this mean?
    753     int nValues = 0;                            // # data points
    754     float rangeMin = 0.0;                       // Exclude data below this
    755     float rangeMax = 0.0;                       // Exclude date above this
    756 
    757     // Determine if the number of data points exceed a threshold which will
    758     // cause to generate robust stats, as opposed to exact stats.
    759     if (myVector->n > stats->sampleLimit) {
    760         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    761         psStats *stats2 = NULL;
    762         // Calculate the robust quartiles.
    763         stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
    764         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    765 
    766         // Store the robust quartiles into the sample quartile members.
    767         stats->sampleUQ = stats2->robustUQ;
    768         stats->sampleLQ = stats2->robustLQ;
    769 
    770         // Free temporary data buffers.
    771         psFree(stats2);
    772 
    773         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    774         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    775 
    776         return;
    777     }
    778 
    779     // Determine how many data points fit inside this min/max range
    780     // and are not maxed, IF the maskVector is not NULL>
    781     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    782 
    783     // Allocate temporary vectors for the data.
    784     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    785     unsortedVector->n = unsortedVector->nalloc;
    786     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
    787     sortedVector->n = sortedVector->nalloc;
    788 
    789     // Determine if we must only use data points within a min/max range.
    790     if (stats->options & PS_STAT_USE_RANGE) {
    791         rangeMin = stats->min;
    792         rangeMax = stats->max;
    793         // Store all non-masked data points within the min/max range
    794         // into the temporary vectors.
    795         count = 0;
    796         if (maskVector != NULL) {
    797             for (i=0;i<myVector->n;i++) {
    798                 if (!(maskVal & maskVector->data.U8[i]) &&
    799                         (rangeMin <= myVector->data.F32[i]) &&
    800                         (myVector->data.F32[i] <= rangeMax)) {
    801                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    802                 }
    803             }
    804         } else {
    805             for (i=0;i<myVector->n;i++) {
    806                 if ((rangeMin <= myVector->data.F32[i]) &&
    807                         (myVector->data.F32[i] <= rangeMax)) {
    808                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    809                 }
    810             }
    811         }
    812     } else {
    813         // Store all non-masked data points into the temporary vectors.
    814         count = 0;
    815         if (maskVector != NULL) {
    816             for (i=0;i<myVector->n;i++) {
    817                 if (!(maskVal & maskVector->data.U8[i])) {
    818                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    819                 }
    820             }
    821         } else {
    822             for (i=0;i<myVector->n;i++) {
    823                 unsortedVector->data.F32[i] = myVector->data.F32[i];
    824             }
    825         }
    826     }
    827 
    828     // Sort the temporary vectors.
    829     psVectorSort(sortedVector, unsortedVector);
    830 
    831     // Calculate the quartile points exactly.
    832     stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
    833     stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    834 
    835     // Free the temporary data structures.
    836     psFree(unsortedVector);
    837     psFree(sortedVector);
    838     // NOTE: This is the
    839 }
    840 
    841 
    842 /******************************************************************************
    843 p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
    844 stdev of the input vector.
    845 Inputs
    846     myVector
    847     maskVector
    848     maskVal
    849     stats
    850 Returns
    851     NULL
    852  
    853 NOTE: the mean is always calculated exactly.  Robust means are never
    854 calculated in this routine.
    855  *****************************************************************************/
    856 void p_psVectorSampleStdev(const psVector *restrict myVector,
    857                            const psVector *restrict maskVector,
    858                            unsigned int maskVal,
    859                            psStats *stats)
    860 {
    861     int i = 0;                                  // Loop index variable
    862     int countInt = 0;                           // # of data points being used
    863     float countFloat = 0.0;                     // # of data points being used
    864     float mean = 0.0;                           // The mean
    865     float diff = 0.0;                           // Used in calculating stdev
    866     float sumSquares = 0.0;                     // temporary variable
    867     float sumDiffs = 0.0;                       // temporary variable
    868     float rangeMin = 0.0;                       // Exclude data below this
    869     float rangeMax = 0.0;                       // Exclude date above this
    870 
    871     // This procedure requires the mean.  If it has not been already
    872     // calculated, then call p_psVectorSampleMean()
    873     if (0 != isnan(stats->sampleMean)) {
    874         p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    875     }
    876     mean = stats->sampleMean;
    877 
    878     if (stats->options & PS_STAT_USE_RANGE) {
    879         if (maskVector != NULL) {
    880             for (i=0;i<myVector->n;i++) {
    881                 if (!(maskVal & maskVector->data.U8[i]) &&
    882                         (rangeMin <= myVector->data.F32[i]) &&
    883                         (myVector->data.F32[i] <= rangeMax)) {
    884                     diff = myVector->data.F32[i] - mean;
    885                     sumSquares+= (diff * diff);
    886                     sumDiffs+= diff;
    887                     countInt++;
    888                 }
    889             }
    890         } else {
    891             for (i=0;i<myVector->n;i++) {
    892                 if ((rangeMin <= myVector->data.F32[i]) &&
    893                         (myVector->data.F32[i] <= rangeMax)) {
    894                     diff = myVector->data.F32[i] - mean;
    895                     sumSquares+= (diff * diff);
    896                     sumDiffs+= diff;
    897                     countInt++;
    898                 }
    899             }
    900             countInt = myVector->n;
    901         }
    902     } else {
    903         if (maskVector != NULL) {
    904             for (i=0;i<myVector->n;i++) {
    905                 if (!(maskVal & maskVector->data.U8[i])) {
    906                     diff = myVector->data.F32[i] - mean;
    907                     sumSquares+= (diff * diff);
    908                     sumDiffs+= diff;
    909                     countInt++;
    910                 }
    911             }
    912         } else {
    913             for (i=0;i<myVector->n;i++) {
    914                 diff = myVector->data.F32[i] - mean;
    915                 sumSquares+= (diff * diff);
    916                 sumDiffs+= diff;
    917                 countInt++;
    918             }
    919             countInt = myVector->n;
    920         }
    921     }
    922     countFloat = (float) countInt;
    923 
    924     #ifdef DARWIN
    925 
    926     stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
    927                                         sumDiffs/countFloat))/ (countFloat-1));
    928     #else
    929 
    930     stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
    931                                  sumDiffs/countFloat))/ (countFloat-1));
    932     #endif
    933 }
    934 
    935 /******************************************************************************
    936 p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
    937 clipped stats (mean or stdev) of the input vector.
    938  
    939 Inputs
    940     myVector
    941     maskVector
    942     maskVal
    943     stats
    944 Returns
    945     NULL
    946  *****************************************************************************/
    947 void p_psVectorClippedStats(const psVector *restrict myVector,
    948                             const psVector *restrict maskVector,
    949                             unsigned int maskVal,
    950                             psStats *stats)
    951 {
    952     int i = 0;                                  // Loop index variable
    953     int j = 0;                                  // Loop index variable
    954     float clippedMean = 0.0;                    // self-explanatory
    955     float clippedStdev = 0.0;                   // self-explanatory
    956     float oldStanMean = 0.0;                    // Temporary variable
    957     float oldStanStdev = 0.0;                   // Temporary variable
    958     psVector *tmpMask = NULL;                   // Temporary vector
    959 
    960     // Endure that stats->clipIter is within the proper range.
    961     if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
    962             (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
    963         psAbort(__func__, "Unallowed value for clipIter (%d).\n",
    964                 stats->clipIter);
    965     }
    966 
    967     // Endure that stats->clipSigma is within the proper range.
    968     if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
    969             (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
    970         psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
    971                 stats->clipSigma);
    972     }
    973 
    974     // We allocate a temporary mask vector since during the iterative
    975     // steps that follow, we will be masking off additional data points.
    976     // However, we do no want to modify the original mask vector.
    977     tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    978     tmpMask->n = myVector->n;
    979 
    980     // If we were called with a mask vector, then initialize the temporary
    981     // mask vector with those values.
    982     if (maskVector != NULL) {
    983         for (i=0;i<tmpMask->n;i++) {
    984             tmpMask->data.U8[i] = maskVector->data.U8[i];
    985         }
    986     }
    987 
    988     // 1. Compute the sample median.
    989     // NOTE: This seems odd.  Verify with IfA that we want to calculate the
    990     // median here, not the mean.
    991     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
    992 
    993     // 2. Compute the sample standard deviation.
    994     p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
    995 
    996     // 3. Use the sample median as the first estimator of the mean X.
    997     clippedMean = stats->sampleMean;
    998 
    999     // 4. Use the sample stdev as the first estimator of the mean stdev.
    1000     clippedStdev = stats->sampleStdev;
    1001 
    1002     // Must save the old sampleMean and sampleStdev since the following code
    1003     // block overwrites them.
    1004     oldStanMean = stats->sampleMean;
    1005     oldStanStdev = stats->sampleStdev;
    1006 
    1007     // 5. Repeat N times:
    1008     for (i=0;i<stats->clipIter;i++) {
    1009         for (j=0;j<myVector->n;j++) {
    1010             // a) Exclude all values x_i for which |x_i - x| > K * stdev
    1011             if (fabs(myVector->data.F32[j] - clippedMean) >
    1012                     (stats->clipSigma * clippedStdev)) {
    1013                 tmpMask->data.U8[i] = 0xff;
    1014             }
    1015             // b) compute new mean and stdev
    1016             p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
    1017             p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
    1018 
    1019             // c) Use the new mean for x
    1020             clippedMean = stats->sampleMean;
    1021 
    1022             // d) Use the new stdev for stdev
    1023             clippedStdev = stats->sampleStdev;
    1024         }
    1025     }
    1026     stats->sampleMean = oldStanMean;
    1027     stats->sampleStdev= oldStanStdev;
    1028 
    1029     // 7. The last calcuated value of x is the cliped mean.
    1030     if (stats->options & PS_STAT_CLIPPED_MEAN) {
    1031         stats->clippedMean = clippedMean;
    1032     }
    1033 
    1034     // 8. The last calcuated value of stdev is the cliped stdev.
    1035     if (stats->options & PS_STAT_CLIPPED_STDEV) {
    1036         stats->clippedStdev = clippedStdev;
    1037     }
    1038 
    1039     psFree(tmpMask);
    1040 }
    1041 
    1042 
    1043 /******************************************************************************
    1044 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
    1045 calculates a variety of robust stat measures:
    1046         PS_STAT_ROBUST_MEAN
    1047         PS_STAT_ROBUST_MEDIAN
    1048         PS_STAT_ROBUST_MODE
    1049         PS_STAT_ROBUST_STDEV
    1050         PS_STAT_ROBUST_QUARTILE
    1051     I have included all that computation in a single function, as opposed to
    1052     breaking it across several functions for one primary reason: they all
    1053     require the same basic initial processing steps (calculate the histogram,
    1054     etc.) however there is no currently defined means, in the SDRS, to
    1055     specify this computation as a separate step.  If the above robust stat
    1056     measures were calcualted in separate functiosn, then much of the initial
    1057     processing would be duplicated.
    1058 Inputs
    1059     myVector
    1060     maskVector
    1061     maskVal
    1062     stats
    1063 Returns
    1064     NULL
    1065  *****************************************************************************/
    1066 void p_psVectorRobustStats(const psVector *restrict myVector,
    1067                            const psVector *restrict maskVector,
    1068                            unsigned int maskVal,
    1069                            psStats *stats)
    1070 {
    1071     psHistogram *robustHistogram = NULL;
    1072     float binSize = 0.0;                        // Size of the histogram bins
    1073     float sigmaE = 0.0;
    1074     int LQBinNum = -1;                          // Bin num for lower quartile
    1075     int UQBinNum = -1;                          // Bin num for upper quartile
    1076     int i = 0;                                  // Loop index variable
    1077     int maxBinNum = 0;
    1078     int maxBinCount = 0;
    1079     float dL = 0.0;
    1080     int numBins = 0;
    1081     psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
    1082 
    1083     // NOTE: The SDRS states that the sample quartiles must be used to
    1084     // determine the initial bin sizes.  However, the sample quartiles are
    1085     // calculated based on a full sort of the data set, regardless of the
    1086     // size of the data set.  We should consult with IfA to ensure that this
    1087     // is really required.
    1088     /*
    1089         if (isnan(stats->sampleUQ) ||
    1090             isnan(stats->sampleLQ)) {
    1091             stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
    1092             p_psVectorSampleQuartiles(myVector,
    1093                                       maskVector,
    1094                                       maskVal,
    1095                                       stats);
    1096         }
    1097     */
    1098     // Compute the initial bin size of the robust histogram.
    1099     p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
    1100     binSize = stats->clippedStdev / 10.0f;
    1101 
    1102     // If stats->clippedStdev == 0.0, then all data elements have the same
    1103     // value.  Therefore, we can set the appropiate results and return.
    1104     if (fabs(binSize) <= FLT_EPSILON) {
    1105         if (stats->options & PS_STAT_ROBUST_MEAN) {
    1106             stats->robustMean = stats->clippedMean;
    1107         }
    1108         if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1109             stats->robustMedian = stats->clippedMean;
    1110         }
    1111         if  (stats->options & PS_STAT_ROBUST_MODE) {
    1112             stats->robustMode = stats->clippedMean;
    1113         }
    1114         if  (stats->options & PS_STAT_ROBUST_STDEV) {
    1115             stats->robustStdev = 0.0;
    1116         }
    1117         if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1118             stats->robustUQ = stats->clippedMean;
    1119             stats->robustLQ = stats->clippedMean;
    1120         }
    1121         psFree(tmpStats);
    1122         psFree(robustHistogram);
    1123         return;
    1124     }
    1125 
    1126     // Detemine minimum and maximum values in the data vector.
    1127     if (isnan(stats->min)) {
    1128         p_psVectorMin(myVector, maskVector, maskVal, stats);
    1129     }
    1130     if (isnan(stats->max)) {
    1131         p_psVectorMax(myVector, maskVector, maskVal, stats);
    1132     }
    1133 
    1134     // Create the histogram structure (yes, 2 is necessary, not 1).  Also,
    1135     // if we get here, we know that binSize != 0.0.
    1136     numBins = 2 + (int) ((stats->max - stats->min) / binSize);
    1137 
    1138     robustHistogram = psHistogramAlloc(stats->min,
    1139                                        stats->max,
    1140                                        numBins);
    1141     // Populate the histogram array.
    1142     psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
    1143 
    1144     // Smooth the histogram.
    1145     p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f);
    1146 
    1147     LQBinNum = -1;
    1148     UQBinNum = -1;
    1149     for (i=0;i<robustHistogram->nums->n;i++) {
    1150         if ((robustHistogram->nums->data.S32[i] <= stats->sampleLQ) &&
    1151                 (stats->sampleLQ <= robustHistogram->nums->data.S32[i])) {
    1152             LQBinNum = i;
    1153         }
    1154 
    1155         if ((robustHistogram->nums->data.S32[i] <= stats->sampleUQ) &&
    1156                 (stats->sampleUQ <= robustHistogram->nums->data.S32[i])) {
    1157             UQBinNum = i;
    1158         }
    1159     }
    1160 
    1161     // Determine the bin with the peak value in the range LQ to UQ.
    1162     maxBinNum = LQBinNum;
    1163     maxBinCount = robustHistogram->nums->data.S32[maxBinNum];
    1164     for (i=LQBinNum;i<=UQBinNum;i++) {
    1165         if (robustHistogram->nums->data.S32[i] > maxBinCount) {
    1166             maxBinNum = i;
    1167             maxBinCount = robustHistogram->nums->data.S32[i];
    1168         }
    1169     }
    1170 
    1171     dL = (stats->robustUQ - stats->robustLQ) / 8.0;
    1172 
    1173     // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    1174     // NOTE: This step is dependent on the functions in psMinimize.c being
    1175     // implemented.  Currently, they are not.
    1176 
    1177     if (stats->options & PS_STAT_ROBUST_MEAN) {
    1178         stats->robustMean = 0.0;
    1179     }
    1180 
    1181     if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1182         stats->robustMedian = 0.0;
    1183     }
    1184     if  (stats->options & PS_STAT_ROBUST_MODE) {
    1185         stats->robustMode = maxBinNum;
    1186     }
    1187     if  (stats->options & PS_STAT_ROBUST_STDEV) {
    1188         stats->robustStdev = 0.0;
    1189     }
    1190     if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1191         stats->robustUQ = 0.0;
    1192         stats->robustLQ = 0.0;
    1193     }
    1194     stats->robustNfit = 0.0;
    1195     stats->robustN50 = 0.0;
    1196     psFree(tmpStats);
    1197     psFree(robustHistogram);
    1198 }
    1199 
    12001758
    12011759/******************************************************************************
    … …  
    12291787    }
    12301788
    1231     // Ensure that the data is of type PS_TYPE_F32.  Eventually, more data
    1232     // types will be implemented.
    1233     if (in->type.type != PS_TYPE_F32) {
    1234         psAbort(__func__,
    1235                 "Only data type PS_TYPE_F32 is currently supported (0x%x).",
    1236                 in->type.type);
    1237     }
    1238 
    1239     // Ensure that the mask vector is of the proper size and type.
     1789    PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
    12401790    if (mask != NULL) {
    1241         if (in->n != mask->n) {
    1242             psAbort(__func__,
    1243                     "Vector data and vector mask are of different sizes.");
    1244         }
    1245         if (mask->type.type != PS_TYPE_U8) {
    1246             psAbort(__func__,
    1247                     "Vector mask must be type PS_TYPE_U8");
    1248         }
     1791        PS_CHECK_NULL_VECTOR(mask);
     1792        PS_CHECK_EMPTY_VECTOR(mask);
     1793        PS_CHECK_VECTOR_SIZE_EQUAL(mask, in);
     1794        PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8);
    12491795    }
    12501796
    … …  
    12811827            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
    12821828        p_psVectorRobustStats(in, mask, maskVal, stats);
     1829
     1830        printf("HMMMM stats->robustMode is %f\n", stats->robustMode);
    12831831    }
    12841832

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

    r1233 r1277  
    55#include <float.h>
    66#include <math.h>
    7 
     7/*****************************************************************************/
     8/* INCLUDE FILES                                                             */
     9/*****************************************************************************/
    810#include "psMemory.h"
     11#include "psImage.h"
    912#include "psVector.h"
    1013#include "psTrace.h"
    … …  
    1215#include "psAbort.h"
    1316#include "psStats.h"
    14 
    15 #include "float.h"
    16 #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
    17 // will use robust statistical methods.
    18 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
    19 #define CLIPPED_NUM_ITER_LB 1
    20 #define CLIPPED_NUM_ITER_UB 10
    21 #define CLIPPED_SIGMA_LB 1.0
    22 #define CLIPPED_SIGMA_UB 10.0
    23 #define true 1
    24 #define false 0
    25 #define MYMAXFLOAT HUGE
    26 
    27 #ifndef DOXYGEN
     17#include "psSort.h"
     18#include "psMinimize.h"
     19#include "psFunctions.h"
     20=======
     21
     22    /*****************************************************************************/
     23    /* DEFINE STATEMENTS                                                         */
     24    /*****************************************************************************/
     25    #define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
     26    // will use robust statistical methods.
     27    #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
     28    #define CLIPPED_NUM_ITER_LB 1
     29    #define CLIPPED_NUM_ITER_UB 10
     30    #define CLIPPED_SIGMA_LB 1.0
     31    #define CLIPPED_SIGMA_UB 10.0
     32    #define true 1
     33    #define false 0
     34    #define MYMAXFLOAT HUGE
     35    #define MAX_ITERATIONS 20
     36
     37    #ifndef DOXYGEN
     38    void p_psVectorRobustStats(const psVector *restrict myVector,
     39                               const psVector *restrict maskVector,
     40                               unsigned int maskVal,
     41                               psStats *stats);
     42#endif
     43
     44/** Preprocessor macro to generate error on an incorrect type */
     45#define PS_CHECK_VECTOR_TYPE(NAME, TYPE)                                                          \
     46if (NAME->type.type != TYPE) {                                                         \
     47    psError(__func__,"Invalid operation: %s has incorrect type.", #NAME);                          \
     48}
     49
     50/** Preprocessor macro to generate error on a NULL vector */
     51#define PS_CHECK_NULL_VECTOR(NAME)                                                          \
     52if (NAME == NULL || NAME->data.V == NULL) {                                                         \
     53    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     54}
     55
     56/** Preprocessor macro to generate error for zero length vector */
     57#define PS_CHECK_EMPTY_VECTOR(NAME)                                                          \
     58if (NAME->n < 1) {                                                                                  \
     59    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                             \
     60}
     61
     62/** Preprocessor macro to generate error on differing size vectors */
     63#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2)                                                          \
     64if (VEC1->n != VEC2->n) {               \
     65    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     66}
     67
     68/*****************************************************************************/
     69/* TYPE DEFINITIONS                                                          */
     70/*****************************************************************************/
     71
     72/*****************************************************************************/
     73/* GLOBAL VARIABLES                                                          */
     74/*****************************************************************************/
     75
     76// None
     77
     78/*****************************************************************************/
     79/* FILE STATIC VARIABLES                                                     */
     80/*****************************************************************************/
     81
     82// None
     83
     84/*****************************************************************************/
     85/* FUNCTION IMPLEMENTATION - LOCAL                                           */
     86/*****************************************************************************/
     87
     88/*****************************************************************************
     89p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
     90function that simply prints a vector to STDOUT.  Used primarily for
     91debugging.
     92 *****************************************************************************/
     93
     94void p_psVectorPrint(psVector *myVector,
     95                     psVector *maskVector,
     96                     unsigned int maskVal,
     97                     psStats *stats)
     98{
     99    int i = 0;                                  // Loop index variable.
     100
     101    for (i=0;i<myVector->n;i++) {
     102        if (maskVector != NULL)
     103            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
     104        else
     105            printf("Element %d is %f\n", i, myVector->data.F32[i]);
     106    }
     107}
     108
     109/******************************************************************************
     110 ******************************************************************************
     111 ******************************************************************************
     112    MISC PRIVATE STATISTICAL FUNCTIONS
     113 
     114    NOTE: it is assumed that any call to these statistical functions will
     115    have been preceded by a call to the psVectorStats() function.  Various
     116    sanity tests will only be performed in psVectorStats().
     117        Is the mask vector the same length as the data vector?
     118        Is the mask vector of type PS_TYPE_U8?
     119        Is the stats data structure NULL?
     120        Is the in data structure NULL?
     121        Is the in data structure of type PS_TYPE_F32?
     122 ******************************************************************************
     123 ******************************************************************************
     124 *****************************************************************************/
     125
     126
     127/******************************************************************************
     128p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
     129mean of the input vector.
     130Inputs
     131Returns
     132    NULL
     133ASSUMPTION: the mean is always calculated exactly.  Robust means are never
     134calculated in this routine.
     135 *****************************************************************************/
     136void p_psVectorSampleMean(const psVector *restrict myVector,
     137                          const psVector *restrict maskVector,
     138                          unsigned int maskVal,
     139                          psStats *stats)
     140{
     141    int i = 0;                                  // Loop index variable
     142    float mean = 0.0;                           // The mean
     143    int count = 0;                              // # of points in this mean?
     144    float rangeMin = 0.0;                       // Exclude data below this
     145    float rangeMax = 0.0;                       // Exclude date above this
     146
     147    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
     148    // loop.
     149    if (stats->options & PS_STAT_USE_RANGE) {
     150        rangeMin = stats->min;
     151        rangeMax = stats->max;
     152        if (maskVector != NULL) {
     153            for (i=0;i<myVector->n;i++) {
     154                // Check if the data is with the specified range
     155                if (!(maskVal & maskVector->data.U8[i]) &&
     156                        (rangeMin <= myVector->data.F32[i]) &&
     157                        (myVector->data.F32[i] <= rangeMax)) {
     158                    mean+= myVector->data.F32[i];
     159                    count++;
     160                }
     161            }
     162            mean/= (float) count;
     163        } else {
     164            for (i=0;i<myVector->n;i++) {
     165                if ((rangeMin <= myVector->data.F32[i]) &&
     166                        (myVector->data.F32[i] <= rangeMax)) {
     167                    mean+= myVector->data.F32[i];
     168                    count++;
     169                }
     170            }
     171            mean/= (float) count;
     172        }
     173    } else {
     174        if (maskVector != NULL) {
     175            for (i=0;i<myVector->n;i++) {
     176                if (!(maskVal & maskVector->data.U8[i])) {
     177                    mean+= myVector->data.F32[i];
     178                    count++;
     179                }
     180            }
     181            mean/= (float) count;
     182        } else {
     183            for (i=0;i<myVector->n;i++) {
     184                mean+= myVector->data.F32[i];
     185            }
     186            mean/= (float) myVector->n;
     187        }
     188    }
     189
     190    stats->sampleMean = mean;
     191}
     192
     193/******************************************************************************
     194p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     195max of the input vector.
     196Inputs
     197    myVector
     198    maskVector
     199    maskVal
     200    stats
     201Returns
     202    NULL
     203 *****************************************************************************/
     204void p_psVectorMax(const psVector *restrict myVector,
     205                   const psVector *restrict maskVector,
     206                   unsigned int maskVal,
     207                   psStats *stats)
     208{
     209    int i = 0;                                  // Loop index variable
     210    float max = -MYMAXFLOAT;                    // The calculated maximum
     211    float rangeMin = 0.0;                       // Exclude data below this
     212    float rangeMax = 0.0;                       // Exclude date above this
     213
     214    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
     215    if (stats->options & PS_STAT_USE_RANGE) {
     216        rangeMin = stats->min;
     217        rangeMax = stats->max;
     218        if (maskVector != NULL) {
     219            for (i=0;i<myVector->n;i++) {
     220                if (!(maskVal & maskVector->data.U8[i])) {
     221                    if ((myVector->data.F32[i] > max) &&
     222                            (rangeMin <= myVector->data.F32[i]) &&
     223                            (myVector->data.F32[i] <= rangeMax)) {
     224                        max = myVector->data.F32[i];
     225                    }
     226                }
     227            }
     228        } else {
     229            for (i=0;i<myVector->n;i++) {
     230                if ((myVector->data.F32[i] > max) &&
     231                        (rangeMin <= myVector->data.F32[i]) &&
     232                        (myVector->data.F32[i] <= rangeMax)) {
     233                    max = myVector->data.F32[i];
     234                }
     235            }
     236        }
     237    } else {
     238        if (maskVector != NULL) {
     239            for (i=0;i<myVector->n;i++) {
     240                if (!(maskVal & maskVector->data.U8[i])) {
     241                    if (myVector->data.F32[i] > max) {
     242                        max = myVector->data.F32[i];
     243                    }
     244                }
     245            }
     246        } else {
     247            for (i=0;i<myVector->n;i++) {
     248                if (myVector->data.F32[i] > max) {
     249                    max = myVector->data.F32[i];
     250                }
     251            }
     252        }
     253    }
     254
     255    stats->max = max;
     256}
     257
     258/******************************************************************************
     259p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
     260minimum of the input vector.
     261Inputs
     262    myVector
     263    maskVector
     264    maskVal
     265    stats
     266Returns
     267    NULL
     268 *****************************************************************************/
     269void p_psVectorMin(const psVector *restrict myVector,
     270                   const psVector *restrict maskVector,
     271                   unsigned int maskVal,
     272                   psStats *stats)
     273{
     274    int i = 0;                                  // Loop index variable
     275    float min = MYMAXFLOAT;                     // The calculated maximum
     276    float rangeMin = 0.0;                       // Exclude data below this
     277    float rangeMax = 0.0;                       // Exclude date above this
     278
     279    if (stats->options & PS_STAT_USE_RANGE) {
     280        rangeMin = stats->min;
     281        rangeMax = stats->max;
     282        if (maskVector != NULL) {
     283            for (i=0;i<myVector->n;i++) {
     284                if (!(maskVal & maskVector->data.U8[i])) {
     285                    if ((myVector->data.F32[i] < min) &&
     286                            (rangeMin <= myVector->data.F32[i]) &&
     287                            (myVector->data.F32[i] <= rangeMax)) {
     288                        min = myVector->data.F32[i];
     289                    }
     290                }
     291            }
     292        } else {
     293            for (i=0;i<myVector->n;i++) {
     294                if ((myVector->data.F32[i] < min) &&
     295                        (rangeMin <= myVector->data.F32[i]) &&
     296                        (myVector->data.F32[i] <= rangeMax)) {
     297                    min = myVector->data.F32[i];
     298                }
     299            }
     300        }
     301    } else {
     302        if (maskVector != NULL) {
     303            for (i=0;i<myVector->n;i++) {
     304                if (!(maskVal & maskVector->data.U8[i])) {
     305                    if (myVector->data.F32[i] < min) {
     306                        min = myVector->data.F32[i];
     307                    }
     308                }
     309            }
     310        } else {
     311            for (i=0;i<myVector->n;i++) {
     312                if (myVector->data.F32[i] < min) {
     313                    min = myVector->data.F32[i];
     314                }
     315            }
     316        }
     317    }
     318
     319    stats->min = min;
     320}
     321
     322/******************************************************************************
     323p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
     324number of non-masked pixels in the vector that fall within the min/max
     325range, if given.
     326Inputs
     327    myVector
     328    maskVector
     329    maskVal
     330    stats
     331Returns
     332    NULL
     333 *****************************************************************************/
     334int p_psVectorNValues(const psVector *restrict myVector,
     335                      const psVector *restrict maskVector,
     336                      unsigned int maskVal,
     337                      psStats *stats)
     338{
     339    int i = 0;                                  // Loop index variable
     340    int numData = 0;                            // The number of data points
     341    float rangeMin = 0.0;                       // Exclude data below this
     342    float rangeMax = 0.0;                       // Exclude date above this
     343
     344    if (stats->options & PS_STAT_USE_RANGE) {
     345        rangeMin = stats->min;
     346        rangeMax = stats->max;
     347        if (maskVector != NULL) {
     348            for (i=0;i<myVector->n;i++) {
     349                if (!(maskVal & maskVector->data.U8[i]) &&
     350                        (rangeMin <= myVector->data.F32[i]) &&
     351                        (myVector->data.F32[i] <= rangeMax)) {
     352                    numData++;
     353                }
     354            }
     355        } else {
     356            for (i=0;i<myVector->n;i++) {
     357                if ((rangeMin <= myVector->data.F32[i]) &&
     358                        (myVector->data.F32[i] <= rangeMax)) {
     359                    numData++;
     360                }
     361            }
     362        }
     363    } else {
     364        rangeMin = stats->min;
     365        rangeMax = stats->max;
     366        if (maskVector != NULL) {
     367            for (i=0;i<myVector->n;i++) {
     368                if (!(maskVal & maskVector->data.U8[i])) {
     369                    numData++;
     370                }
     371            }
     372        } else {
     373            numData = myVector->n;
     374        }
     375    }
     376    return(numData);
     377}
     378
     379
     380
     381/******************************************************************************
     382p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
     383median of the input vector.
     384Inputs
     385    myVector
     386    maskVector
     387    maskVal
     388    stats
     389Returns
     390    NULL
     391 *****************************************************************************/
     392void p_psVectorSampleMedian(const psVector *restrict myVector,
     393                            const psVector *restrict maskVector,
     394                            unsigned int maskVal,
     395                            psStats *stats)
     396{
     397    psVector *unsortedVector = NULL;            // Temporary vector
     398    psVector *sortedVector = NULL;              // Temporary vector
     399    int i = 0;                                  // Loop index variable
     400    int count = 0;                              // # of points in this mean?
     401    int nValues = 0;                            // # of points in vector
     402    float rangeMin = 0.0;                       // Exclude data below this
     403    float rangeMax = 0.0;                       // Exclude date above this
     404    psStats *stats2 = NULL;                     // Temporary stats structure
     405
     406
     407    // Determine if the number of data points exceed a threshold which will
     408    // cause to generate robust stats, as opposed to exact stats.
     409
     410    if (myVector->n > stats->sampleLimit) {
     411        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     412
     413        // Calculate the robust quartiles.
     414        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     415        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     416
     417        // Store the robust quartiles into the sample quartile members.
     418        stats->sampleMedian = stats2->robustMedian;
     419
     420        // Free temporary data buffers.
     421        psFree(stats2);
     422
     423        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     424        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     425
     426        return;
     427    }
     428
     429    // Determine how many data points fit inside this min/max range
     430    // and are not masked, IF the maskVector is not NULL>
     431    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
     432
     433    // Allocate temporary vectors for the data.
     434    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     435    unsortedVector->n = unsortedVector->nalloc;
     436
     437    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     438    sortedVector->n = sortedVector->nalloc;
     439
     440    // Determine if we must only use data points within a min/max range.
     441    if (stats->options & PS_STAT_USE_RANGE) {
     442        rangeMin = stats->min;
     443        rangeMax = stats->max;
     444
     445        // Store all non-masked data points within the min/max range
     446        // into the temporary vectors.
     447        count = 0;
     448        if (maskVector != NULL) {
     449            for (i=0;i<myVector->n;i++) {
     450                if (!(maskVal & maskVector->data.U8[i]) &&
     451                        (rangeMin <= myVector->data.F32[i]) &&
     452                        (myVector->data.F32[i] <= rangeMax)) {
     453                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
     454                }
     455            }
     456        } else {
     457            for (i=0;i<myVector->n;i++) {
     458                if ((rangeMin <= myVector->data.F32[i]) &&
     459                        (myVector->data.F32[i] <= rangeMax)) {
     460                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     461                }
     462            }
     463        }
     464    } else {
     465        // Store all non-masked data points into the temporary vectors.
     466        count = 0;
     467        if (maskVector != NULL) {
     468            for (i=0;i<myVector->n;i++) {
     469                if (!(maskVal & maskVector->data.U8[i])) {
     470                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     471                }
     472            }
     473        } else {
     474            for (i=0;i<myVector->n;i++) {
     475                unsortedVector->data.F32[i] = myVector->data.F32[i];
     476            }
     477        }
     478    }
     479    // Sort the temporary vectors.
     480    psVectorSort(sortedVector, unsortedVector);
     481
     482    // Calculate the median exactly.
     483    if (0 == (nValues % 2)) {
     484        stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
     485                                     sortedVector->data.F32[nValues/2]);
     486    } else {
     487        stats->sampleMedian = sortedVector->data.F32[nValues/2];
     488    }
     489
     490    // Free the temporary data structures.
     491    psFree(unsortedVector);
     492    psFree(sortedVector);
     493}
     494
     495/******************************************************************************
     496    This routine smoothes the data in the input robustHistogram with a
     497    Gaussian of width sigma.
     498 *****************************************************************************/
     499psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
     500                                       float sigma)
     501{
     502    int i = 0;                                  // Loop index variable
     503    int j = 0;                                  // Loop index variable
     504    float denom = 0.0;                           // Temporary variable
     505    float expo = 0.0;                            // Temporary variable
     506    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
     507    psVector *smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
     508
     509    for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
     510        if (fabs(sigma) >= FLT_EPSILON) {
     511            // If sigma does not equal zero, then we use Gaussian smoothing.
     512            #ifdef  DARWIN
     513            denom = (float) sqrt(2.0 * M_PI * sigma * sigma);
     514            #else
     515
     516            denom = sqrtf(2.0 * M_PI * sigma * sigma);
     517            #endif
     518
     519            expo = - (float) ((i-GAUSS_WIDTH) * (i-GAUSS_WIDTH));
     520            expo/= (2.0 * sigma * sigma);
     521            gaussianCoefs[i] = exp(expo/denom);
     522
     523            // NOTE: Gaussian smoothing just isn't working with low sigma
     524            // values.  The problem is that the Gaussian coefficients are
     525            // all zero, except for the middle coefficient, which is exactly
     526            // one.  Therefore, I'm using boxcar smoothing.
     527            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float) GAUSS_WIDTH));
     528            //            printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
     529        } else {
     530            /* If sigma equals zero (all pixels have the same value)
     531             * the above code will divide by zero.  Therefore, we don't need
     532             * to smooth the data.
     533             */
     534            for(i=0;i<robustHistogram->nums->n;i++) {
     535                smooth->data.F32[i] = (float) robustHistogram->nums->data.S32[i];
     536            }
     537            return(smooth);
     538        }
     539    }
     540
     541    for(i=0;i<robustHistogram->nums->n;i++) {
     542        smooth->data.F32[i] = 0.0;
     543        for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
     544            if (((j+i) >= 0) && ((j+i) < smooth->n)) {
     545                smooth->data.F32[i]+= (gaussianCoefs[j+GAUSS_WIDTH] *
     546                                       (float) robustHistogram->nums->data.S32[j+i]);
     547            }
     548        }
     549    }
     550    return(smooth);
     551}
     552
     553/******************************************************************************
     554p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
     555the upper and/or lower quartiles of the input vector.
     556Inputs
     557    myVector
     558    maskVector
     559    maskVal
     560    stats
     561Returns
     562    NULL
     563 *****************************************************************************/
     564void p_psVectorSampleQuartiles(const psVector *restrict myVector,
     565                               const psVector *restrict maskVector,
     566                               unsigned int maskVal,
     567                               psStats *stats)
     568{
     569    psVector *unsortedVector = NULL;            // Temporary vector
     570    psVector *sortedVector = NULL;              // Temporary vector
     571    int i = 0;                                  // Loop index variable
     572    int count = 0;                              // # of points in this mean?
     573    int nValues = 0;                            // # data points
     574    float rangeMin = 0.0;                       // Exclude data below this
     575    float rangeMax = 0.0;                       // Exclude date above this
     576
     577    // Determine if the number of data points exceed a threshold which will
     578    // cause to generate robust stats, as opposed to exact stats.
     579    if (myVector->n > stats->sampleLimit) {
     580        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     581        psStats *stats2 = NULL;
     582        // Calculate the robust quartiles.
     583        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
     584        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     585
     586        // Store the robust quartiles into the sample quartile members.
     587        stats->sampleUQ = stats2->robustUQ;
     588        stats->sampleLQ = stats2->robustLQ;
     589
     590        // Free temporary data buffers.
     591        psFree(stats2);
     592
     593        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     594        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     595
     596        return;
     597    }
     598
     599    // Determine how many data points fit inside this min/max range
     600    // and are not maxed, IF the maskVector is not NULL>
     601    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
     602
     603    // Allocate temporary vectors for the data.
     604    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
     605    unsortedVector->n = unsortedVector->nalloc;
     606    sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
     607    sortedVector->n = sortedVector->nalloc;
     608
     609    // Determine if we must only use data points within a min/max range.
     610    if (stats->options & PS_STAT_USE_RANGE) {
     611        rangeMin = stats->min;
     612        rangeMax = stats->max;
     613        // Store all non-masked data points within the min/max range
     614        // into the temporary vectors.
     615        count = 0;
     616        if (maskVector != NULL) {
     617            for (i=0;i<myVector->n;i++) {
     618                if (!(maskVal & maskVector->data.U8[i]) &&
     619                        (rangeMin <= myVector->data.F32[i]) &&
     620                        (myVector->data.F32[i] <= rangeMax)) {
     621                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     622                }
     623            }
     624        } else {
     625            for (i=0;i<myVector->n;i++) {
     626                if ((rangeMin <= myVector->data.F32[i]) &&
     627                        (myVector->data.F32[i] <= rangeMax)) {
     628                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     629                }
     630            }
     631        }
     632    } else {
     633        // Store all non-masked data points into the temporary vectors.
     634        count = 0;
     635        if (maskVector != NULL) {
     636            for (i=0;i<myVector->n;i++) {
     637                if (!(maskVal & maskVector->data.U8[i])) {
     638                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
     639                }
     640            }
     641        } else {
     642            for (i=0;i<myVector->n;i++) {
     643                unsortedVector->data.F32[i] = myVector->data.F32[i];
     644            }
     645        }
     646    }
     647
     648    // Sort the temporary vectors.
     649    psVectorSort(sortedVector, unsortedVector);
     650
     651    // Calculate the quartile points exactly.
     652    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     653    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
     654
     655    // Free the temporary data structures.
     656    psFree(unsortedVector);
     657    psFree(sortedVector);
     658    // NOTE: This is the
     659}
     660
     661
     662/******************************************************************************
     663p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
     664stdev of the input vector.
     665Inputs
     666    myVector
     667    maskVector
     668    maskVal
     669    stats
     670Returns
     671    NULL
     672 
     673NOTE: the mean is always calculated exactly.  Robust means are never
     674calculated in this routine.
     675 *****************************************************************************/
     676void p_psVectorSampleStdev(const psVector *restrict myVector,
     677                           const psVector *restrict maskVector,
     678                           unsigned int maskVal,
     679                           psStats *stats)
     680{
     681    int i = 0;                                  // Loop index variable
     682    int countInt = 0;                           // # of data points being used
     683    float countFloat = 0.0;                     // # of data points being used
     684    float mean = 0.0;                           // The mean
     685    float diff = 0.0;                           // Used in calculating stdev
     686    float sumSquares = 0.0;                     // temporary variable
     687    float sumDiffs = 0.0;                       // temporary variable
     688    float rangeMin = 0.0;                       // Exclude data below this
     689    float rangeMax = 0.0;                       // Exclude date above this
     690
     691    // This procedure requires the mean.  If it has not been already
     692    // calculated, then call p_psVectorSampleMean()
     693    if (0 != isnan(stats->sampleMean)) {
     694        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     695    }
     696    mean = stats->sampleMean;
     697
     698    if (stats->options & PS_STAT_USE_RANGE) {
     699        if (maskVector != NULL) {
     700            for (i=0;i<myVector->n;i++) {
     701                if (!(maskVal & maskVector->data.U8[i]) &&
     702                        (rangeMin <= myVector->data.F32[i]) &&
     703                        (myVector->data.F32[i] <= rangeMax)) {
     704                    diff = myVector->data.F32[i] - mean;
     705                    sumSquares+= (diff * diff);
     706                    sumDiffs+= diff;
     707                    countInt++;
     708                }
     709            }
     710        } else {
     711            for (i=0;i<myVector->n;i++) {
     712                if ((rangeMin <= myVector->data.F32[i]) &&
     713                        (myVector->data.F32[i] <= rangeMax)) {
     714                    diff = myVector->data.F32[i] - mean;
     715                    sumSquares+= (diff * diff);
     716                    sumDiffs+= diff;
     717                    countInt++;
     718                }
     719            }
     720            countInt = myVector->n;
     721        }
     722    } else {
     723        if (maskVector != NULL) {
     724            for (i=0;i<myVector->n;i++) {
     725                if (!(maskVal & maskVector->data.U8[i])) {
     726                    diff = myVector->data.F32[i] - mean;
     727                    sumSquares+= (diff * diff);
     728                    sumDiffs+= diff;
     729                    countInt++;
     730                }
     731            }
     732        } else {
     733            for (i=0;i<myVector->n;i++) {
     734                diff = myVector->data.F32[i] - mean;
     735                sumSquares+= (diff * diff);
     736                sumDiffs+= diff;
     737                countInt++;
     738            }
     739            countInt = myVector->n;
     740        }
     741    }
     742    countFloat = (float) countInt;
     743
     744    #ifdef DARWIN
     745
     746    stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
     747                                        sumDiffs/countFloat))/ (countFloat-1));
     748    #else
     749
     750    stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
     751                                 sumDiffs/countFloat))/ (countFloat-1));
     752    #endif
     753}
     754
     755/******************************************************************************
     756p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
     757clipped stats (mean or stdev) of the input vector.
     758 
     759Inputs
     760    myVector
     761    maskVector
     762    maskVal
     763    stats
     764Returns
     765    NULL
     766 *****************************************************************************/
     767void p_psVectorClippedStats(const psVector *restrict myVector,
     768                            const psVector *restrict maskVector,
     769                            unsigned int maskVal,
     770                            psStats *stats)
     771{
     772    int i = 0;                                  // Loop index variable
     773    int j = 0;                                  // Loop index variable
     774    float clippedMean = 0.0;                    // self-explanatory
     775    float clippedStdev = 0.0;                   // self-explanatory
     776    float oldStanMean = 0.0;                    // Temporary variable
     777    float oldStanStdev = 0.0;                   // Temporary variable
     778    psVector *tmpMask = NULL;                   // Temporary vector
     779
     780    // Endure that stats->clipIter is within the proper range.
     781    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
     782            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
     783        psAbort(__func__, "Unallowed value for clipIter (%d).\n",
     784                stats->clipIter);
     785    }
     786
     787    // Endure that stats->clipSigma is within the proper range.
     788    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
     789            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
     790        psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
     791                stats->clipSigma);
     792    }
     793
     794    // We allocate a temporary mask vector since during the iterative
     795    // steps that follow, we will be masking off additional data points.
     796    // However, we do no want to modify the original mask vector.
     797    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
     798    tmpMask->n = myVector->n;
     799
     800    // If we were called with a mask vector, then initialize the temporary
     801    // mask vector with those values.
     802    if (maskVector != NULL) {
     803        for (i=0;i<tmpMask->n;i++) {
     804            tmpMask->data.U8[i] = maskVector->data.U8[i];
     805        }
     806    }
     807
     808    // 1. Compute the sample median.
     809    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
     810    // median here, not the mean.
     811    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
     812
     813    // 2. Compute the sample standard deviation.
     814    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
     815
     816    // 3. Use the sample median as the first estimator of the mean X.
     817    clippedMean = stats->sampleMean;
     818
     819    // 4. Use the sample stdev as the first estimator of the mean stdev.
     820    clippedStdev = stats->sampleStdev;
     821
     822    // Must save the old sampleMean and sampleStdev since the following code
     823    // block overwrites them.
     824    oldStanMean = stats->sampleMean;
     825    oldStanStdev = stats->sampleStdev;
     826
     827    // 5. Repeat N times:
     828    for (i=0;i<stats->clipIter;i++) {
     829        for (j=0;j<myVector->n;j++) {
     830            // a) Exclude all values x_i for which |x_i - x| > K * stdev
     831            if (fabs(myVector->data.F32[j] - clippedMean) >
     832                    (stats->clipSigma * clippedStdev)) {
     833                tmpMask->data.U8[i] = 0xff;
     834            }
     835            // b) compute new mean and stdev
     836            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
     837            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
     838
     839            // c) Use the new mean for x
     840            clippedMean = stats->sampleMean;
     841
     842            // d) Use the new stdev for stdev
     843            clippedStdev = stats->sampleStdev;
     844        }
     845
     846    }
     847    stats->sampleMean = oldStanMean;
     848    stats->sampleStdev= oldStanStdev;
     849
     850    // 7. The last calcuated value of x is the cliped mean.
     851    if (stats->options & PS_STAT_CLIPPED_MEAN) {
     852        stats->clippedMean = clippedMean;
     853    }
     854
     855    // 8. The last calcuated value of stdev is the cliped stdev.
     856    if (stats->options & PS_STAT_CLIPPED_STDEV) {
     857        stats->clippedStdev = clippedStdev;
     858    }
     859
     860    psFree(tmpMask);
     861}
     862
     863/*****************************************************************************
     864 *****************************************************************************/
     865void p_psNormalizeVector(psVector *myData)
     866{
     867    float min = (float) HUGE;
     868    float max = (float) -HUGE;
     869    float range = 0.0;
     870    int i = 0;
     871
     872    for (i=0;i<myData->n;i++) {
     873        if (myData->data.F32[i] < min) {
     874            min = myData->data.F32[i];
     875        }
     876        if (myData->data.F32[i] > max) {
     877            max = myData->data.F32[i];
     878        }
     879    }
     880
     881    range = max - min;
     882    for (i=0;i<myData->n;i++) {
     883        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
     884    }
     885}
     886
     887
     888float p_psGaussian(const psVector *restrict myData,
     889                   const psVector *restrict myParams)
     890{
     891    float x = myData->data.F32[0];
     892    float mean = myParams->data.F32[0];
     893    float stdev = myParams->data.F32[1];
     894    float tmp = exp(-((x-mean) * (x-mean)) / (2.0 * stdev * stdev));
     895    tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
     896
     897    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     898    return(tmp);
     899}
     900
     901float p_psGaussianDeriv(const psVector *restrict myData,
     902                        const psVector *restrict myParams,
     903                        int whichParam)
     904{
     905    float x = myData->data.F32[0];
     906    float mean = myParams->data.F32[0];
     907    float stdev = myParams->data.F32[1];
     908    float tmp = 0.0;
     909
     910    if (whichParam == 0) {
     911        // Return the derivative w.r.t. the mean.
     912        tmp = (x - mean) * p_psGaussian(myData, myParams);
     913        tmp/= (stdev * stdev);
     914    } else if (whichParam == 1) {
     915        // Return the derivative w.r.t. the stdev.
     916        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
     917        tmp/= (stdev * stdev * stdev);
     918    }
     919    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
     920
     921    return(tmp);
     922}
     923
     924
     925float p_psQuadratic(const psVector *restrict myParams,
     926                    const psVector *restrict myCoords)
     927{
     928    float x = myCoords->data.F32[0];
     929    float A = myParams->data.F32[0];
     930    float B = myParams->data.F32[1];
     931    float C = myParams->data.F32[2];
     932    float tmp = 0.0;
     933
     934    tmp = (A * x * x) + (B * x) + C;
     935    return(tmp);
     936}
     937
     938float p_psQuadraticDeriv(const psVector *restrict myParams,
     939                         const psVector *restrict myCoords,
     940                         int whichParamDeriv)
     941{
     942    float x = myCoords->data.F32[0];
     943    float tmp = 0.0;
     944
     945    if (whichParamDeriv == 0) {
     946        tmp = x * x;
     947    } else if (whichParamDeriv == 1) {
     948        tmp = x;
     949    } else if (whichParamDeriv == 2) {
     950        tmp = 1.0;
     951    }
     952
     953    return(tmp);
     954}
     955
     956/******************************************************************************
     957p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
     958as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
     959polynomials here) and a range of x-values for which it is defined:
     960[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
     961that f(x) == midpoint.  This functions uses a binary-search algorithm on the
     962range and assumes that the polnomial is monotonically increasing within that
     963range.
     964 *****************************************************************************/
     965float p_ps1DPolyMedian(psPolynomial1D *myPoly,
     966                       float rangeLow,
     967                       float rangeHigh,
     968                       float getThisValue)
     969{
     970    int numIterations=0;
     971    float midpoint = 0.0;
     972    float oldMidpoint = 1.0;
     973    float f = 0.0;
     974
     975    //  printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
     976
     977    while (numIterations < MAX_ITERATIONS) {
     978        midpoint = (rangeHigh + rangeLow) / 2.0;
     979        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
     980            return(midpoint);
     981        }
     982        oldMidpoint = midpoint;
     983
     984        f = psEvalPolynomial1D(midpoint, myPoly);
     985        //        printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
     986        if (fabs(f - getThisValue) <= FLT_EPSILON) {
     987            return(midpoint);
     988        }
     989
     990        if (f > getThisValue) {
     991            rangeHigh = midpoint;
     992        } else {
     993            rangeLow = midpoint;
     994        }
     995        numIterations++;
     996    }
     997    return(midpoint);
     998}
     999
     1000/******************************************************************************
     1001p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
     1002as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
     1003polynomials here) and a range of x-values for which it is defined:
     1004[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
     1005that f(x) == midpoint.  This functions uses a binary-search algorithm on the
     1006range and assumes that the polnomial is monotonically increasing within that
     1007range.
     1008 *****************************************************************************/
     1009float p_psFitQuadratic(psHistogram *robustHistogram,
     1010                       int binNum,
     1011                       float fitFloat)
     1012{
     1013    /*
     1014        if ((binNum > 0) &&
     1015            (binNum < (robustHistogram->nums->n+1))) {
     1016            x->data.F64[0] = (double) 0.5 *
     1017                      (robustHistogram->bounds->data.F32[binNum-1] +
     1018                       robustHistogram->bounds->data.F32[binNum]);
     1019            x->data.F64[1] = (double) 0.5 *
     1020                      (robustHistogram->bounds->data.F32[binNum] +
     1021                       robustHistogram->bounds->data.F32[binNum+1]);
     1022            x->data.F64[2] = (double) 0.5 *
     1023                      (robustHistogram->bounds->data.F32[binNum+1] +
     1024                       robustHistogram->bounds->data.F32[binNum+2]);
     1025     
     1026            y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1];
     1027            y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum];
     1028            y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1];
     1029     
     1030            if (!((y->data.F64[0] <= fitFloat) &&
     1031                 (fitFloat <= y->data.F64[2]))) {
     1032                psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
     1033            }
     1034     
     1035            yErr->data.F64[0] = 1.0;
     1036            yErr->data.F64[1] = 1.0;
     1037            yErr->data.F64[2] = 1.0;
     1038     
     1039            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1040            return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat);
     1041         } else {
     1042            return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] +
     1043                          robustHistogram->bounds->data.F32[binNum]));
     1044        }
     1045    */
     1046    return(0.0);
     1047}
     1048
     1049#define PS_PRINT_VECTOR(NAME) \
     1050for (int my_i=0;my_i<NAME->n;my_i++) { \
     1051    printf("%s->data.F32[%d] is %f\n", #NAME, my_i, NAME->data.F32[my_i]); \
     1052} \
     1053printf("\n"); \
     1054
     1055/******************************************************************************
     1056p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
     1057calculates a variety of robust stat measures:
     1058PS_STAT_ROBUST_MEAN
     1059PS_STAT_ROBUST_MEDIAN
     1060PS_STAT_ROBUST_MODE
     1061PS_STAT_ROBUST_STDEV
     1062PS_STAT_ROBUST_QUARTILE
     1063I have included all that computation in a single function, as opposed to
     1064breaking it across several functions for one primary reason:
     1065they all
     1066require the same basic initial processing steps (calculate the histogram,
     1067etc.)
     1068however there is no currently defined means, in the SDRS, to
     1069specify this computation as a separate step.  If the above robust stat
     1070measures were calcualted in separate functiosn, then much of the initial
     1071processing would be duplicated.
     1072Inputs
     1073myVector
     1074maskVector
     1075maskVal
     1076stats
     1077Returns
     1078NULL
     1079*****************************************************************************/
    281080void p_psVectorRobustStats(const psVector *restrict myVector,
    291081                           const psVector *restrict maskVector,
    301082                           unsigned int maskVal,
    31                            psStats *stats);
    32 #endif
     1083                           psStats *stats)
     1084{
     1085    psHistogram *robustHistogram = NULL;
     1086    psVector *robustHistogramVector = NULL;
     1087    float binSize = 0.0;                        // Size of the histogram bins
     1088    int LQBinNum = -1;                          // Bin num for lower quartile
     1089    int UQBinNum = -1;                          // Bin num for upper quartile
     1090    int i = 0;                                  // Loop index variable
     1091    int maxBinNum = 0;
     1092    float maxBinCount = 0.0;
     1093    float dL = 0.0;
     1094    int numBins = 0;
     1095    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
     1096    //    psImage *domain;
     1097    //    psVector *errors;
     1098    //    psVector *data;
     1099    //    psVector *initialGuess;
     1100    //    psVector *theParams;
     1101    //    float chiSq=0.0;
     1102    //    float max = -HUGE;
     1103    //    float max_pos;
     1104    float myMean = 0.0;
     1105    float myStdev = 0.0;
     1106    float countFloat = 0.0;
     1107    float diff = 0.0;
     1108    float sumSquares = 0.0;
     1109    float sumDiffs = 0.0;
     1110    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1111    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1112    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1113    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
     1114    psVector *cumulativeRobustSumsFullRange = NULL;
     1115    psVector *cumulativeRobustSumsDlRange = NULL;
     1116    float sumRobust = 0.0;
     1117    float sumN50 = 0.0;
     1118    float sumNfit = 0.0;
     1119    float cumulativeMedian = 0.0;
     1120
     1121    // Compute the initial bin size of the robust histogram.  This is done
     1122    // by computing the clipped standard deviation of the vector, and dividing
     1123    // that by 10.0;
     1124    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
     1125    binSize = tmpStats->clippedStdev / 10.0f;
     1126
     1127    // If stats->clippedStdev == 0.0, then all data elements have the same
     1128    // value.  Therefore, we can set the appropiate results and return.
     1129    if (fabs(binSize) <= FLT_EPSILON) {
     1130        if (stats->options & PS_STAT_ROBUST_MEAN) {
     1131            stats->robustMean = stats->clippedMean;
     1132        }
     1133        if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
     1134            stats->robustMedian = stats->clippedMean;
     1135        }
     1136        if  (stats->options & PS_STAT_ROBUST_MODE) {
     1137            stats->robustMode = stats->clippedMean;
     1138        }
     1139        if  (stats->options & PS_STAT_ROBUST_STDEV) {
     1140            stats->robustStdev = 0.0;
     1141        }
     1142        if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
     1143            stats->robustUQ = stats->clippedMean;
     1144            stats->robustLQ = stats->clippedMean;
     1145        }
     1146        // XXX: Set these to the number of unmasked data points?
     1147        stats->robustNfit = 0.0;
     1148        stats->robustN50 = 0.0;
     1149        psFree(tmpStats);
     1150        return;
     1151    }
     1152
     1153    // Determine minimum and maximum values in the data vector.
     1154    if (isnan(stats->min)) {
     1155        p_psVectorMin(myVector, maskVector, maskVal, stats);
     1156    }
     1157    if (isnan(stats->max)) {
     1158        p_psVectorMax(myVector, maskVector, maskVal, stats);
     1159    }
     1160
     1161    // Create the histogram structure.  NOTE: we can not specify the bin size
     1162    // precisely since the argument to psHistogramAlloc() is the number of
     1163    // bins, not the binSize.  Also, if we get here, we know that
     1164    // binSize != 0.0.
     1165    numBins = (int) ((stats->max - stats->min) / binSize);
     1166    robustHistogram = psHistogramAlloc(stats->min,
     1167                                       stats->max,
     1168                                       numBins);
     1169
     1170    // Populate the histogram array.
     1171    psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
     1172
     1173    // Smooth the histogram.
     1174    // XXX: is that the right stdev?
     1175    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
     1176                            tmpStats->clippedStdev/4.0f);
     1177
     1178    // The following was necessary to fit a gaussian to the data, since
     1179    // gaussian functions produce data between 0.0 and 1.0.
     1180    // p_psNormalizeVector(robustHistogramVector);
     1181
     1182    /**************************************************************************
     1183    Determine the lower/upper quartiles.
     1184    **************************************************************************/
     1185    // We define a vector called "cumulativeRobustSums..." where the value at
     1186    // index position i is equal to the sum of bins 0:i.  This will be used
     1187    // now and later in determining the lower/upper quartiles.
     1188    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1189    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
     1190    for (i=1;i<robustHistogramVector->n;i++) {
     1191        cumulativeRobustSumsFullRange->data.F32[i] =
     1192            cumulativeRobustSumsFullRange->data.F32[i-1] +
     1193            robustHistogramVector->data.F32[i];
     1194    }
     1195    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1];
     1196
     1197    // Determine the bin number containing the lower quartile point.
     1198    LQBinNum = -1;
     1199    for (i=0;i<cumulativeRobustSumsFullRange->n;i++) {
     1200        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust/4.0)) {
     1201            LQBinNum = i;
     1202            break;
     1203        }
     1204    }
     1205
     1206    // Determine the bin number containing the upper quartile point.
     1207    UQBinNum = -1;
     1208    for (i=cumulativeRobustSumsFullRange->n-1;i>=0;i--) {
     1209        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0*sumRobust/4.0)) {
     1210            UQBinNum = i;
     1211            break;
     1212        }
     1213    }
     1214
     1215    if ((LQBinNum == -1) ||
     1216            (UQBinNum == -1)) {
     1217        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
     1218    }
     1219    /**************************************************************************
     1220    Determine the mode in the range LQ:UQ.
     1221    **************************************************************************/
     1222    // Determine the bin with the peak value in the range LQ to UQ.
     1223    maxBinNum = LQBinNum;
     1224    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
     1225    sumN50 = (float) robustHistogram->nums->data.S32[LQBinNum];
     1226    for (i=LQBinNum+1;i<=UQBinNum;i++) {
     1227        if (robustHistogramVector->data.F32[i] > maxBinCount) {
     1228            maxBinNum = i;
     1229            maxBinCount = robustHistogramVector->data.F32[i];
     1230        }
     1231        sumN50+= (float) robustHistogram->nums->data.S32[i];
     1232    }
     1233    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
     1234    dL = (UQBinNum - LQBinNum) / 4;
     1235
     1236    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
     1237
     1238    /**************************************************************************
     1239    Determine the mean/stdev for the bins in the range mode-dL to mode+dL
     1240    **************************************************************************/
     1241    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1242    for (i=0;i<robustHistogramVector->n;i++) {
     1243        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
     1244    }
     1245    sumNfit = 0.0;
     1246    cumulativeMedian = 0.0;
     1247    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1248        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1249            cumulativeRobustSumsDlRange->data.F32[i] =
     1250                cumulativeRobustSumsDlRange->data.F32[i-1] +
     1251                robustHistogramVector->data.F32[i];
     1252            cumulativeMedian+= robustHistogramVector->data.F32[i];
     1253            sumNfit+= (float) robustHistogramVector->data.S32[i];
     1254        }
     1255    }
     1256
     1257    // Calculate the mean of the smoothed robust histogram in the range
     1258    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     1259    // that bin (this is a non-exact approximation).
     1260    myMean = 0.0;
     1261    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1262        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1263            myMean+= (robustHistogramVector->data.F32[i]) * 0.5 *
     1264                     (robustHistogram->bounds->data.F32[i+1] +
     1265                      robustHistogram->bounds->data.F32[i]);
     1266            countFloat+= robustHistogramVector->data.F32[i];
     1267        }
     1268    }
     1269    myMean/= countFloat;
     1270
     1271    // Calculate the stdev of the smoothed robust histogram in the range
     1272    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     1273    // that bin.
     1274    for (i=maxBinNum-dL;i<=maxBinNum+dL;i++) {
     1275        if ((0 <= i) && (i < robustHistogramVector->n)) {
     1276            diff = (0.5 * (robustHistogram->bounds->data.F32[i+1] +
     1277                           robustHistogram->bounds->data.F32[i])) - myMean;
     1278            sumSquares+= diff * diff * robustHistogramVector->data.F32[i];
     1279            sumDiffs+= diff * robustHistogramVector->data.F32[i];
     1280        }
     1281    }
     1282    myStdev = sqrt((sumSquares-(sumDiffs * sumDiffs/countFloat))/ (countFloat-1));
     1283
     1284    /**************************************************************************
     1285    Set the appropriate members in the output stats struct.
     1286    **************************************************************************/
     1287    if (stats->options & PS_STAT_ROBUST_MEAN) {
     1288        stats->robustMean = myMean;
     1289    }
     1290
     1291    if  (stats->options & PS_STAT_ROBUST_MODE) {
     1292        stats->robustMode = 0.5 *
     1293                            (robustHistogram->bounds->data.F32[maxBinNum] +
     1294                             robustHistogram->bounds->data.F32[maxBinNum+1]);
     1295    }
     1296
     1297    if  (stats->options & PS_STAT_ROBUST_STDEV) {
     1298        stats->robustStdev = myStdev;
     1299    }
     1300
     1301    // To determine the median (and later, the lower/upper quartile), we fit
     1302    // a quadratic to the three bins surrounding the bin containing the median.
     1303    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
     1304    // the cumulative number of data points in all bins up to, and including,
     1305    // this bin.  We then solve the quadratic for
     1306
     1307    if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
     1308        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n-1))) {
     1309            x->data.F64[0] = (double) 0.5 *
     1310                             (robustHistogram->bounds->data.F32[maxBinNum-1] +
     1311                              robustHistogram->bounds->data.F32[maxBinNum]);
     1312            x->data.F64[1] = (double) 0.5 *
     1313                             (robustHistogram->bounds->data.F32[maxBinNum] +
     1314                              robustHistogram->bounds->data.F32[maxBinNum+1]);
     1315            x->data.F64[2] = (double) 0.5 *
     1316                             (robustHistogram->bounds->data.F32[maxBinNum+1] +
     1317                              robustHistogram->bounds->data.F32[maxBinNum+2]);
     1318
     1319            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum-1];
     1320            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
     1321            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum+1];
     1322
     1323            // Ensure that cumulativeMedian/2 is actually within the range of the bins
     1324            // we are using.
     1325            cumulativeMedian*= 0.5;
     1326            if (!((y->data.F64[0] <= cumulativeMedian) &&
     1327                    (cumulativeMedian <= y->data.F64[2]))) {
     1328                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
     1329                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
     1330            }
     1331            // XXX: yErr is not currently used by psGetArrayPolynomial().  We
     1332            // may have to set this meaningfully later.
     1333            yErr->data.F64[0] = 1.0;
     1334            yErr->data.F64[1] = 1.0;
     1335            yErr->data.F64[2] = 1.0;
     1336
     1337            // Determine the coefficients of the polynomial.
     1338            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1339            // Call p_ps1DPolyMedian(), which does a binary search on the
     1340            // polynomial, looking for the value x such that
     1341            // f(x) = cumulativeMedian.
     1342            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0],
     1343                                                   x->data.F64[2], cumulativeMedian);
     1344        } else {
     1345            // If the mode is the first/last histogram bin, then simply use
     1346            // the midpoint of that bin.
     1347            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum+1] +
     1348                                         robustHistogram->bounds->data.F32[maxBinNum]);
     1349        }
     1350    }
     1351
     1352    // The lower/upper quartile calculations are very similar to the median
     1353    // calculations.  We fit a quadratic to the array containing the
     1354    // cumulative data points in each bin, then search for x such that
     1355    // f(x) equals the lower/upper quartile exactly.
     1356    //
     1357    if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
     1358        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n-1];
     1359
     1360        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n-1))) {
     1361            x->data.F64[0] = (double) 0.5 *
     1362                             (robustHistogram->bounds->data.F32[LQBinNum-1] +
     1363                              robustHistogram->bounds->data.F32[LQBinNum]);
     1364            x->data.F64[1] = (double) 0.5 *
     1365                             (robustHistogram->bounds->data.F32[LQBinNum] +
     1366                              robustHistogram->bounds->data.F32[LQBinNum+1]);
     1367            x->data.F64[2] = (double) 0.5 *
     1368                             (robustHistogram->bounds->data.F32[LQBinNum+1] +
     1369                              robustHistogram->bounds->data.F32[LQBinNum+2]);
     1370
     1371            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum-1];
     1372            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
     1373            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum+1];
     1374
     1375            if (!((y->data.F64[0] <= (countFloat/4.0)) &&
     1376                    ((countFloat/4.0) <= y->data.F64[2]))) {
     1377                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
     1378            }
     1379
     1380            yErr->data.F64[0] = 1.0;
     1381            yErr->data.F64[1] = 1.0;
     1382            yErr->data.F64[2] = 1.0;
     1383
     1384            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1385            stats->robustLQ = p_ps1DPolyMedian(myPoly,
     1386                                               x->data.F64[0],
     1387                                               x->data.F64[2],
     1388                                               countFloat/4.0);
     1389
     1390        } else {
     1391            // If the LQ is the first/last histogram bin, then simply use
     1392            // the midpoint of that bin.
     1393            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum+1] +
     1394                                     robustHistogram->bounds->data.F32[LQBinNum]);
     1395        }
     1396
     1397        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n-1))) {
     1398            x->data.F64[0] = (double) 0.5 *
     1399                             (robustHistogram->bounds->data.F32[UQBinNum-1] +
     1400                              robustHistogram->bounds->data.F32[UQBinNum]);
     1401            x->data.F64[1] = (double) 0.5 *
     1402                             (robustHistogram->bounds->data.F32[UQBinNum] +
     1403                              robustHistogram->bounds->data.F32[UQBinNum+1]);
     1404            x->data.F64[2] = (double) 0.5 *
     1405                             (robustHistogram->bounds->data.F32[UQBinNum+1] +
     1406                              robustHistogram->bounds->data.F32[UQBinNum+2]);
     1407
     1408            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum-1];
     1409            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
     1410            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum+1];
     1411
     1412            if (!((y->data.F64[0] <= (3.0 * countFloat/4.0)) &&
     1413                    ((3.0 * countFloat/4.0) <= y->data.F64[2]))) {
     1414                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
     1415            }
     1416
     1417            yErr->data.F64[0] = 1.0;
     1418            yErr->data.F64[1] = 1.0;
     1419            yErr->data.F64[2] = 1.0;
     1420
     1421            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1422            stats->robustUQ = p_ps1DPolyMedian(myPoly,
     1423                                               x->data.F64[0],
     1424                                               x->data.F64[2],
     1425                                               3.0*countFloat/4.0);
     1426        } else {
     1427            // If the UQ is the first/last histogram bin, then simply use
     1428            // the midpoint of that bin.
     1429            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum+1] +
     1430                                     robustHistogram->bounds->data.F32[UQBinNum]);
     1431        }
     1432    }
     1433    stats->robustNfit = sumNfit;
     1434    stats->robustN50 = sumN50;
     1435
     1436    psFree(x);
     1437    psFree(y);
     1438    psFree(yErr);
     1439    psFree(tmpStats);
     1440    psFree(robustHistogram);
     1441    psFree(myPoly);
     1442    psFree(cumulativeRobustSumsFullRange);
     1443    psFree(cumulativeRobustSumsDlRange);
     1444}
     1445
     1446
     1447
     1448/*
     1449void p_ps_FitTheGaussian()
     1450{
     1451 
     1452    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     1453    // NOTE: This code uses the psMinimize.c functions to perform
     1454    // the fit.  It doesn't quite work, 100% of the time.  For the time being
     1455    // I am commenting this code out, and replacing it with code which
     1456    // calculates the mean directly on the robustHistogram.
     1457 
     1458    domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
     1459    errors = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1460    data = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1461    initialGuess = psVectorAlloc(2, PS_TYPE_F32);
     1462 
     1463    max = -HUGE;
     1464    max_pos = -1;
     1465    for (i=0;i<robustHistogramVector->n;i++) {
     1466        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
     1467                                        robustHistogram->bounds->data.F32[i]);
     1468        data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
     1469        errors->data.F32[i] = 1.0;
     1470        //printf("DATA (%.2f, %.2f)\n", domain->data.F32[i][0], data->data.F32[i]);
     1471        if (data->data.F32[i] > max) {
     1472            max = data->data.F32[i];
     1473            max_pos = domain->data.F32[i][0];
     1474        }
     1475    }
     1476 
     1477    initialGuess->data.F32[0] = max_pos;
     1478    initialGuess->data.F32[1] = 1.0;
     1479    printf("Initial (mean. stdev) is (%.3f, %.3f)\n", initialGuess->data.F32[0], initialGuess->data.F32[1]);
     1480 
     1481    printf("Calling psMinimizeChi2()\n");
     1482    theParams = psMinimizeChi2(p_psGaussian,
     1483                               p_psGaussianDeriv,
     1484                               domain,
     1485                               data,
     1486                               errors,
     1487                               initialGuess,
     1488                               NULL,
     1489                               &chiSq);
     1490    printf("Called psMinimizeChi2()\n");
     1491    printf("p_psVectorRobustStats() is (%f, %f)\n", theParams->data.F32[0], theParams->data.F32[1]);
     1492}
     1493*/
     1494
     1495
     1496/*****************************************************************************/
     1497/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     1498/*****************************************************************************/
    331499
    341500static void histogramFree(psHistogram *myHist);
    … …  
    921558        return(NULL);
    931559    }
     1560    if (n < 0) {
     1561        psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n);
     1562    }
    941563
    951564    // NOTE: Verify that this is the correct action.
    … …  
    1071576    // Calculate the bounds for each bin.
    1081577    binSize = (upper - lower) / (float) n;
     1578    // NOTE: Is the following necessary?  It prevents the max data point
     1579    // from being in a non-existant bin.
     1580    binSize+= FLT_EPSILON;
    1091581    for (i=0;i<n+1;i++) {
    1101582        newHist->bounds->data.F32[i] = lower + (binSize * (float) i);
    … …  
    1271599
    1281600/******************************************************************************
    129 psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure
    130 with the specifed bounds.  The number of elements in the bounds vector is n.
    131 Therefore, the number of bins is n-1.
     1601psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure
     1602with the specifed bounds.
     1603 
    1321604Input:
    1331605    bounds
    … …  
    2381710
    2391711    numBins = out->nums->n;
    240 
    2411712    for (i=0;i<in->n;i++) {
    2421713        // Check if this pixel is masked, and if so, skip it.
    … …  
    2561727                if (out->uniform == true) {
    2571728                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    258 
    2591729                    binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    2601730                                    binSize);
     1731
     1732                    // NOTE: This next if-statement really shouldn't be necessary.
     1733                    // However, do to numerical lack of precision, we occasionally
     1734                    // produce a binNum outside the range of bins.
     1735                    if (binNum >= out->nums->n) {
     1736                        binNum = out->nums->n-1;
     1737                    }
     1738
    2611739                    (out->nums->data.S32[binNum])++;
     1740
    2621741                    // If this is a non-uniform histogram, determining the correct
    2631742                    // bin number requires a bit more work.
    … …  
    2771756    return(out);
    2781757}
    279 
    280 /*****************************************************************************
    281 p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
    282 function that simply prints a vector to STDOUT.  Used primarily for
    283 debugging.
    284  *****************************************************************************/
    285 
    286 void p_psVectorPrint(psVector *myVector,
    287                      psVector *maskVector,
    288                      unsigned int maskVal,
    289                      psStats *stats)
    290 {
    291     int i = 0;                                  // Loop index variable.
    292 
    293     for (i=0;i<myVector->n;i++) {
    294         if (maskVector != NULL)
    295             printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
    296         else
    297             printf("Element %d is %f\n", i, myVector->data.F32[i]);
    298     }
    299 }
    300 
    301 
    302 /******************************************************************************
    303  ******************************************************************************
    304  ******************************************************************************
    305     MISC PRIVATE STATISTICAL FUNCTIONS
    306  
    307     NOTE: it is assumed that any call to these statistical functions will
    308     have been preceded by a call to the psVectorStats() function.  Various
    309     sanity tests will only be performed in psVectorStats().
    310         Is the mask vector the same length as the data vector?
    311         Is the mask vector of type PS_TYPE_U8?
    312         Is the stats data structure NULL?
    313         Is the in data structure NULL?
    314         Is the in data structure of type PS_TYPE_F32?
    315  ******************************************************************************
    316  ******************************************************************************
    317  *****************************************************************************/
    318 
    319 
    320 /******************************************************************************
    321 p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
    322 mean of the input vector.
    323 Inputs
    324 Returns
    325     NULL
    326 ASSUMPTION: the mean is always calculated exactly.  Robust means are never
    327 calculated in this routine.
    328  *****************************************************************************/
    329 void p_psVectorSampleMean(const psVector *restrict myVector,
    330                           const psVector *restrict maskVector,
    331                           unsigned int maskVal,
    332                           psStats *stats)
    333 {
    334     int i = 0;                                  // Loop index variable
    335     float mean = 0.0;                           // The mean
    336     int count = 0;                              // # of points in this mean?
    337     float rangeMin = 0.0;                       // Exclude data below this
    338     float rangeMax = 0.0;                       // Exclude date above this
    339 
    340     // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
    341     // loop.
    342     if (stats->options & PS_STAT_USE_RANGE) {
    343         rangeMin = stats->min;
    344         rangeMax = stats->max;
    345         if (maskVector != NULL) {
    346             for (i=0;i<myVector->n;i++) {
    347                 // Check if the data is with the specified range
    348                 if (!(maskVal & maskVector->data.U8[i]) &&
    349                         (rangeMin <= myVector->data.F32[i]) &&
    350                         (myVector->data.F32[i] <= rangeMax)) {
    351                     mean+= myVector->data.F32[i];
    352                     count++;
    353                 }
    354             }
    355             mean/= (float) count;
    356         } else {
    357             for (i=0;i<myVector->n;i++) {
    358                 if ((rangeMin <= myVector->data.F32[i]) &&
    359                         (myVector->data.F32[i] <= rangeMax)) {
    360                     mean+= myVector->data.F32[i];
    361                     count++;
    362                 }
    363             }
    364             mean/= (float) count;
    365         }
    366     } else {
    367         if (maskVector != NULL) {
    368             for (i=0;i<myVector->n;i++) {
    369                 if (!(maskVal & maskVector->data.U8[i])) {
    370                     mean+= myVector->data.F32[i];
    371                     count++;
    372                 }
    373             }
    374             mean/= (float) count;
    375         } else {
    376             for (i=0;i<myVector->n;i++) {
    377                 mean+= myVector->data.F32[i];
    378             }
    379             mean/= (float) myVector->n;
    380         }
    381     }
    382 
    383     stats->sampleMean = mean;
    384 }
    385 
    386 /******************************************************************************
    387 p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
    388 max of the input vector.
    389 Inputs
    390     myVector
    391     maskVector
    392     maskVal
    393     stats
    394 Returns
    395     NULL
    396  *****************************************************************************/
    397 void p_psVectorMax(const psVector *restrict myVector,
    398                    const psVector *restrict maskVector,
    399                    unsigned int maskVal,
    400                    psStats *stats)
    401 {
    402     int i = 0;                                  // Loop index variable
    403     float max = -MYMAXFLOAT;                    // The calculated maximum
    404     float rangeMin = 0.0;                       // Exclude data below this
    405     float rangeMax = 0.0;                       // Exclude date above this
    406 
    407     // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    408     if (stats->options & PS_STAT_USE_RANGE) {
    409         rangeMin = stats->min;
    410         rangeMax = stats->max;
    411         if (maskVector != NULL) {
    412             for (i=0;i<myVector->n;i++) {
    413                 if (!(maskVal & maskVector->data.U8[i])) {
    414                     if ((myVector->data.F32[i] > max) &&
    415                             (rangeMin <= myVector->data.F32[i]) &&
    416                             (myVector->data.F32[i] <= rangeMax)) {
    417                         max = myVector->data.F32[i];
    418                     }
    419                 }
    420             }
    421         } else {
    422             for (i=0;i<myVector->n;i++) {
    423                 if ((myVector->data.F32[i] > max) &&
    424                         (rangeMin <= myVector->data.F32[i]) &&
    425                         (myVector->data.F32[i] <= rangeMax)) {
    426                     max = myVector->data.F32[i];
    427                 }
    428             }
    429         }
    430     } else {
    431         if (maskVector != NULL) {
    432             for (i=0;i<myVector->n;i++) {
    433                 if (!(maskVal & maskVector->data.U8[i])) {
    434                     if (myVector->data.F32[i] > max) {
    435                         max = myVector->data.F32[i];
    436                     }
    437                 }
    438             }
    439         } else {
    440             for (i=0;i<myVector->n;i++) {
    441                 if (myVector->data.F32[i] > max) {
    442                     max = myVector->data.F32[i];
    443                 }
    444             }
    445         }
    446     }
    447 
    448     stats->max = max;
    449 }
    450 
    451 /******************************************************************************
    452 p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
    453 minimum of the input vector.
    454 Inputs
    455     myVector
    456     maskVector
    457     maskVal
    458     stats
    459 Returns
    460     NULL
    461  *****************************************************************************/
    462 void p_psVectorMin(const psVector *restrict myVector,
    463                    const psVector *restrict maskVector,
    464                    unsigned int maskVal,
    465                    psStats *stats)
    466 {
    467     int i = 0;                                  // Loop index variable
    468     float min = MYMAXFLOAT;                     // The calculated maximum
    469     float rangeMin = 0.0;                       // Exclude data below this
    470     float rangeMax = 0.0;                       // Exclude date above this
    471 
    472     if (stats->options & PS_STAT_USE_RANGE) {
    473         rangeMin = stats->min;
    474         rangeMax = stats->max;
    475         if (maskVector != NULL) {
    476             for (i=0;i<myVector->n;i++) {
    477                 if (!(maskVal & maskVector->data.U8[i])) {
    478                     if ((myVector->data.F32[i] < min) &&
    479                             (rangeMin <= myVector->data.F32[i]) &&
    480                             (myVector->data.F32[i] <= rangeMax)) {
    481                         min = myVector->data.F32[i];
    482                     }
    483                 }
    484             }
    485         } else {
    486             for (i=0;i<myVector->n;i++) {
    487                 if ((myVector->data.F32[i] < min) &&
    488                         (rangeMin <= myVector->data.F32[i]) &&
    489                         (myVector->data.F32[i] <= rangeMax)) {
    490                     min = myVector->data.F32[i];
    491                 }
    492             }
    493         }
    494     } else {
    495         if (maskVector != NULL) {
    496             for (i=0;i<myVector->n;i++) {
    497                 if (!(maskVal & maskVector->data.U8[i])) {
    498                     if (myVector->data.F32[i] < min) {
    499                         min = myVector->data.F32[i];
    500                     }
    501                 }
    502             }
    503         } else {
    504             for (i=0;i<myVector->n;i++) {
    505                 if (myVector->data.F32[i] < min) {
    506                     min = myVector->data.F32[i];
    507                 }
    508             }
    509         }
    510     }
    511 
    512     stats->min = min;
    513 }
    514 
    515 /******************************************************************************
    516 p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
    517 number of non-masked pixels in the vector that fall within the min/max
    518 range, if given.
    519 Inputs
    520     myVector
    521     maskVector
    522     maskVal
    523     stats
    524 Returns
    525     NULL
    526  *****************************************************************************/
    527 int p_psVectorNValues(const psVector *restrict myVector,
    528                       const psVector *restrict maskVector,
    529                       unsigned int maskVal,
    530                       psStats *stats)
    531 {
    532     int i = 0;                                  // Loop index variable
    533     int numData = 0;                            // The number of data points
    534     float rangeMin = 0.0;                       // Exclude data below this
    535     float rangeMax = 0.0;                       // Exclude date above this
    536 
    537     if (stats->options & PS_STAT_USE_RANGE) {
    538         rangeMin = stats->min;
    539         rangeMax = stats->max;
    540         if (maskVector != NULL) {
    541             for (i=0;i<myVector->n;i++) {
    542                 if (!(maskVal & maskVector->data.U8[i]) &&
    543                         (rangeMin <= myVector->data.F32[i]) &&
    544                         (myVector->data.F32[i] <= rangeMax)) {
    545                     numData++;
    546                 }
    547             }
    548         } else {
    549             for (i=0;i<myVector->n;i++) {
    550                 if ((rangeMin <= myVector->data.F32[i]) &&
    551                         (myVector->data.F32[i] <= rangeMax)) {
    552                     numData++;
    553                 }
    554             }
    555         }
    556     } else {
    557         rangeMin = stats->min;
    558         rangeMax = stats->max;
    559         if (maskVector != NULL) {
    560             for (i=0;i<myVector->n;i++) {
    561                 if (!(maskVal & maskVector->data.U8[i])) {
    562                     numData++;
    563                 }
    564             }
    565         } else {
    566             numData = myVector->n;
    567         }
    568     }
    569     return(numData);
    570 }
    571 
    572 
    573 
    574 /******************************************************************************
    575 p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
    576 median of the input vector.
    577 Inputs
    578     myVector
    579     maskVector
    580     maskVal
    581     stats
    582 Returns
    583     NULL
    584  *****************************************************************************/
    585 void p_psVectorSampleMedian(const psVector *restrict myVector,
    586                             const psVector *restrict maskVector,
    587                             unsigned int maskVal,
    588                             psStats *stats)
    589 {
    590     psVector *unsortedVector = NULL;            // Temporary vector
    591     psVector *sortedVector = NULL;              // Temporary vector
    592     int i = 0;                                  // Loop index variable
    593     int count = 0;                              // # of points in this mean?
    594     int nValues = 0;                            // # of points in vector
    595     float rangeMin = 0.0;                       // Exclude data below this
    596     float rangeMax = 0.0;                       // Exclude date above this
    597     psStats *stats2 = NULL;                     // Temporary stats structure
    598 
    599 
    600     // Determine if the number of data points exceed a threshold which will
    601     // cause to generate robust stats, as opposed to exact stats.
    602 
    603     if (myVector->n > stats->sampleLimit) {
    604         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    605 
    606         // Calculate the robust quartiles.
    607         stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    608         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    609 
    610         // Store the robust quartiles into the sample quartile members.
    611         stats->sampleMedian = stats2->robustMedian;
    612 
    613         // Free temporary data buffers.
    614         psFree(stats2);
    615 
    616         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    617         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    618 
    619         return;
    620     }
    621 
    622     // Determine how many data points fit inside this min/max range
    623     // and are not masked, IF the maskVector is not NULL>
    624     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    625 
    626     // Allocate temporary vectors for the data.
    627     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    628     unsortedVector->n = unsortedVector->nalloc;
    629 
    630     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
    631     sortedVector->n = sortedVector->nalloc;
    632 
    633     // Determine if we must only use data points within a min/max range.
    634     if (stats->options & PS_STAT_USE_RANGE) {
    635         rangeMin = stats->min;
    636         rangeMax = stats->max;
    637 
    638         // Store all non-masked data points within the min/max range
    639         // into the temporary vectors.
    640         count = 0;
    641         if (maskVector != NULL) {
    642             for (i=0;i<myVector->n;i++) {
    643                 if (!(maskVal & maskVector->data.U8[i]) &&
    644                         (rangeMin <= myVector->data.F32[i]) &&
    645                         (myVector->data.F32[i] <= rangeMax)) {
    646                     unsortedVector->data.F32[count++] = maskVector->data.F32[i];
    647                 }
    648             }
    649         } else {
    650             for (i=0;i<myVector->n;i++) {
    651                 if ((rangeMin <= myVector->data.F32[i]) &&
    652                         (myVector->data.F32[i] <= rangeMax)) {
    653                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    654                 }
    655             }
    656         }
    657     } else {
    658         // Store all non-masked data points into the temporary vectors.
    659         count = 0;
    660         if (maskVector != NULL) {
    661             for (i=0;i<myVector->n;i++) {
    662                 if (!(maskVal & maskVector->data.U8[i])) {
    663                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    664                 }
    665             }
    666         } else {
    667             for (i=0;i<myVector->n;i++) {
    668                 unsortedVector->data.F32[i] = myVector->data.F32[i];
    669             }
    670         }
    671     }
    672     // Sort the temporary vectors.
    673     psVectorSort(sortedVector, unsortedVector);
    674 
    675     // Calculate the median exactly.
    676     if (0 == (nValues % 2)) {
    677         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
    678                                      sortedVector->data.F32[nValues/2]);
    679     } else {
    680         stats->sampleMedian = sortedVector->data.F32[nValues/2];
    681     }
    682 
    683     // Free the temporary data structures.
    684     psFree(unsortedVector);
    685     psFree(sortedVector);
    686 }
    687 
    688 /******************************************************************************
    689     This routine smoothes the data in the input robustHistogram with a
    690     Gaussian of width sigma.
    691  *****************************************************************************/
    692 void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
    693                                   float sigma)
    694 {
    695     int i = 0;                                  // Loop index variable
    696     int j = 0;                                  // Loop index variable
    697     float tmpf = 0.0;                           // Temporary variable
    698     float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
    699 
    700     for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
    701         if (fabs(sigma) <= FLT_EPSILON) {
    702             // If sigma does not equal zero, then we use Gaussian smoothing.
    703             #ifdef  DARWIN
    704             tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma);
    705             #else
    706 
    707             tmpf = sqrtf(2.0f * M_PI * sigma * sigma);
    708             #endif
    709 
    710             gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) *
    711                                              ((float) (i-GAUSS_WIDTH))) /
    712                                             (2.0f * sigma * sigma)) / tmpf;
    713         } else {
    714             /* If sigma equals zero (all pixels have the same value)
    715              * the above code will divide by zero.  Therefore, we don't need
    716              * to smooth the data.
    717              */
    718             return;
    719         }
    720     }
    721 
    722     for(i=0;i<robustHistogram->nums->n;i++) {
    723         for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
    724             if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) {
    725                 robustHistogram->nums->data.S32[j+i]+=
    726                     (gaussianCoefs[j+GAUSS_WIDTH] *
    727                      (float) robustHistogram->nums->data.S32[j+i]);
    728             }
    729         }
    730     }
    731 }
    732 
    733 /******************************************************************************
    734 p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
    735 the upper and/or lower quartiles of the input vector.
    736 Inputs
    737     myVector
    738     maskVector
    739     maskVal
    740     stats
    741 Returns
    742     NULL
    743  *****************************************************************************/
    744 void p_psVectorSampleQuartiles(const psVector *restrict myVector,
    745                                const psVector *restrict maskVector,
    746                                unsigned int maskVal,
    747                                psStats *stats)
    748 {
    749     psVector *unsortedVector = NULL;            // Temporary vector
    750     psVector *sortedVector = NULL;              // Temporary vector
    751     int i = 0;                                  // Loop index variable
    752     int count = 0;                              // # of points in this mean?
    753     int nValues = 0;                            // # data points
    754     float rangeMin = 0.0;                       // Exclude data below this
    755     float rangeMax = 0.0;                       // Exclude date above this
    756 
    757     // Determine if the number of data points exceed a threshold which will
    758     // cause to generate robust stats, as opposed to exact stats.
    759     if (myVector->n > stats->sampleLimit) {
    760         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    761         psStats *stats2 = NULL;
    762         // Calculate the robust quartiles.
    763         stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
    764         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    765 
    766         // Store the robust quartiles into the sample quartile members.
    767         stats->sampleUQ = stats2->robustUQ;
    768         stats->sampleLQ = stats2->robustLQ;
    769 
    770         // Free temporary data buffers.
    771         psFree(stats2);
    772 
    773         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    774         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    775 
    776         return;
    777     }
    778 
    779     // Determine how many data points fit inside this min/max range
    780     // and are not maxed, IF the maskVector is not NULL>
    781     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    782 
    783     // Allocate temporary vectors for the data.
    784     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    785     unsortedVector->n = unsortedVector->nalloc;
    786     sortedVector   = psVectorAlloc(nValues, PS_TYPE_F32);
    787     sortedVector->n = sortedVector->nalloc;
    788 
    789     // Determine if we must only use data points within a min/max range.
    790     if (stats->options & PS_STAT_USE_RANGE) {
    791         rangeMin = stats->min;
    792         rangeMax = stats->max;
    793         // Store all non-masked data points within the min/max range
    794         // into the temporary vectors.
    795         count = 0;
    796         if (maskVector != NULL) {
    797             for (i=0;i<myVector->n;i++) {
    798                 if (!(maskVal & maskVector->data.U8[i]) &&
    799                         (rangeMin <= myVector->data.F32[i]) &&
    800                         (myVector->data.F32[i] <= rangeMax)) {
    801                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    802                 }
    803             }
    804         } else {
    805             for (i=0;i<myVector->n;i++) {
    806                 if ((rangeMin <= myVector->data.F32[i]) &&
    807                         (myVector->data.F32[i] <= rangeMax)) {
    808                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    809                 }
    810             }
    811         }
    812     } else {
    813         // Store all non-masked data points into the temporary vectors.
    814         count = 0;
    815         if (maskVector != NULL) {
    816             for (i=0;i<myVector->n;i++) {
    817                 if (!(maskVal & maskVector->data.U8[i])) {
    818                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
    819                 }
    820             }
    821         } else {
    822             for (i=0;i<myVector->n;i++) {
    823                 unsortedVector->data.F32[i] = myVector->data.F32[i];
    824             }
    825         }
    826     }
    827 
    828     // Sort the temporary vectors.
    829     psVectorSort(sortedVector, unsortedVector);
    830 
    831     // Calculate the quartile points exactly.
    832     stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
    833     stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    834 
    835     // Free the temporary data structures.
    836     psFree(unsortedVector);
    837     psFree(sortedVector);
    838     // NOTE: This is the
    839 }
    840 
    841 
    842 /******************************************************************************
    843 p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
    844 stdev of the input vector.
    845 Inputs
    846     myVector
    847     maskVector
    848     maskVal
    849     stats
    850 Returns
    851     NULL
    852  
    853 NOTE: the mean is always calculated exactly.  Robust means are never
    854 calculated in this routine.
    855  *****************************************************************************/
    856 void p_psVectorSampleStdev(const psVector *restrict myVector,
    857                            const psVector *restrict maskVector,
    858                            unsigned int maskVal,
    859                            psStats *stats)
    860 {
    861     int i = 0;                                  // Loop index variable
    862     int countInt = 0;                           // # of data points being used
    863     float countFloat = 0.0;                     // # of data points being used
    864     float mean = 0.0;                           // The mean
    865     float diff = 0.0;                           // Used in calculating stdev
    866     float sumSquares = 0.0;                     // temporary variable
    867     float sumDiffs = 0.0;                       // temporary variable
    868     float rangeMin = 0.0;                       // Exclude data below this
    869     float rangeMax = 0.0;                       // Exclude date above this
    870 
    871     // This procedure requires the mean.  If it has not been already
    872     // calculated, then call p_psVectorSampleMean()
    873     if (0 != isnan(stats->sampleMean)) {
    874         p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    875     }
    876     mean = stats->sampleMean;
    877 
    878     if (stats->options & PS_STAT_USE_RANGE) {
    879         if (maskVector != NULL) {
    880             for (i=0;i<myVector->n;i++) {
    881                 if (!(maskVal & maskVector->data.U8[i]) &&
    882                         (rangeMin <= myVector->data.F32[i]) &&
    883                         (myVector->data.F32[i] <= rangeMax)) {
    884                     diff = myVector->data.F32[i] - mean;
    885                     sumSquares+= (diff * diff);
    886                     sumDiffs+= diff;
    887                     countInt++;
    888                 }
    889             }
    890         } else {
    891             for (i=0;i<myVector->n;i++) {
    892                 if ((rangeMin <= myVector->data.F32[i]) &&
    893                         (myVector->data.F32[i] <= rangeMax)) {
    894                     diff = myVector->data.F32[i] - mean;
    895                     sumSquares+= (diff * diff);
    896                     sumDiffs+= diff;
    897                     countInt++;
    898                 }
    899             }
    900             countInt = myVector->n;
    901         }
    902     } else {
    903         if (maskVector != NULL) {
    904             for (i=0;i<myVector->n;i++) {
    905                 if (!(maskVal & maskVector->data.U8[i])) {
    906                     diff = myVector->data.F32[i] - mean;
    907                     sumSquares+= (diff * diff);
    908                     sumDiffs+= diff;
    909                     countInt++;
    910                 }
    911             }
    912         } else {
    913             for (i=0;i<myVector->n;i++) {
    914                 diff = myVector->data.F32[i] - mean;
    915                 sumSquares+= (diff * diff);
    916                 sumDiffs+= diff;
    917                 countInt++;
    918             }
    919             countInt = myVector->n;
    920         }
    921     }
    922     countFloat = (float) countInt;
    923 
    924     #ifdef DARWIN
    925 
    926     stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
    927                                         sumDiffs/countFloat))/ (countFloat-1));
    928     #else
    929 
    930     stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
    931                                  sumDiffs/countFloat))/ (countFloat-1));
    932     #endif
    933 }
    934 
    935 /******************************************************************************
    936 p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
    937 clipped stats (mean or stdev) of the input vector.
    938  
    939 Inputs
    940     myVector
    941     maskVector
    942     maskVal
    943     stats
    944 Returns
    945     NULL
    946  *****************************************************************************/
    947 void p_psVectorClippedStats(const psVector *restrict myVector,
    948                             const psVector *restrict maskVector,
    949                             unsigned int maskVal,
    950                             psStats *stats)
    951 {
    952     int i = 0;                                  // Loop index variable
    953     int j = 0;                                  // Loop index variable
    954     float clippedMean = 0.0;                    // self-explanatory
    955     float clippedStdev = 0.0;                   // self-explanatory
    956     float oldStanMean = 0.0;                    // Temporary variable
    957     float oldStanStdev = 0.0;                   // Temporary variable
    958     psVector *tmpMask = NULL;                   // Temporary vector
    959 
    960     // Endure that stats->clipIter is within the proper range.
    961     if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
    962             (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
    963         psAbort(__func__, "Unallowed value for clipIter (%d).\n",
    964                 stats->clipIter);
    965     }
    966 
    967     // Endure that stats->clipSigma is within the proper range.
    968     if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
    969             (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
    970         psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
    971                 stats->clipSigma);
    972     }
    973 
    974     // We allocate a temporary mask vector since during the iterative
    975     // steps that follow, we will be masking off additional data points.
    976     // However, we do no want to modify the original mask vector.
    977     tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    978     tmpMask->n = myVector->n;
    979 
    980     // If we were called with a mask vector, then initialize the temporary
    981     // mask vector with those values.
    982     if (maskVector != NULL) {
    983         for (i=0;i<tmpMask->n;i++) {
    984             tmpMask->data.U8[i] = maskVector->data.U8[i];
    985         }
    986     }
    987 
    988     // 1. Compute the sample median.
    989     // NOTE: This seems odd.  Verify with IfA that we want to calculate the
    990     // median here, not the mean.
    991     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
    992 
    993     // 2. Compute the sample standard deviation.
    994     p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
    995 
    996     // 3. Use the sample median as the first estimator of the mean X.
    997     clippedMean = stats->sampleMean;
    998 
    999     // 4. Use the sample stdev as the first estimator of the mean stdev.
    1000     clippedStdev = stats->sampleStdev;
    1001 
    1002     // Must save the old sampleMean and sampleStdev since the following code
    1003     // block overwrites them.
    1004     oldStanMean = stats->sampleMean;
    1005     oldStanStdev = stats->sampleStdev;
    1006 
    1007     // 5. Repeat N times:
    1008     for (i=0;i<stats->clipIter;i++) {
    1009         for (j=0;j<myVector->n;j++) {
    1010             // a) Exclude all values x_i for which |x_i - x| > K * stdev
    1011             if (fabs(myVector->data.F32[j] - clippedMean) >
    1012                     (stats->clipSigma * clippedStdev)) {
    1013                 tmpMask->data.U8[i] = 0xff;
    1014             }
    1015             // b) compute new mean and stdev
    1016             p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
    1017             p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
    1018 
    1019             // c) Use the new mean for x
    1020             clippedMean = stats->sampleMean;
    1021 
    1022             // d) Use the new stdev for stdev
    1023             clippedStdev = stats->sampleStdev;
    1024         }
    1025     }
    1026     stats->sampleMean = oldStanMean;
    1027     stats->sampleStdev= oldStanStdev;
    1028 
    1029     // 7. The last calcuated value of x is the cliped mean.
    1030     if (stats->options & PS_STAT_CLIPPED_MEAN) {
    1031         stats->clippedMean = clippedMean;
    1032     }
    1033 
    1034     // 8. The last calcuated value of stdev is the cliped stdev.
    1035     if (stats->options & PS_STAT_CLIPPED_STDEV) {
    1036         stats->clippedStdev = clippedStdev;
    1037     }
    1038 
    1039     psFree(tmpMask);
    1040 }
    1041 
    1042 
    1043 /******************************************************************************
    1044 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
    1045 calculates a variety of robust stat measures:
    1046         PS_STAT_ROBUST_MEAN
    1047         PS_STAT_ROBUST_MEDIAN
    1048         PS_STAT_ROBUST_MODE
    1049         PS_STAT_ROBUST_STDEV
    1050         PS_STAT_ROBUST_QUARTILE
    1051     I have included all that computation in a single function, as opposed to
    1052     breaking it across several functions for one primary reason: they all
    1053     require the same basic initial processing steps (calculate the histogram,
    1054     etc.) however there is no currently defined means, in the SDRS, to
    1055     specify this computation as a separate step.  If the above robust stat
    1056     measures were calcualted in separate functiosn, then much of the initial
    1057     processing would be duplicated.
    1058 Inputs
    1059     myVector
    1060     maskVector
    1061     maskVal
    1062     stats
    1063 Returns
    1064     NULL
    1065  *****************************************************************************/
    1066 void p_psVectorRobustStats(const psVector *restrict myVector,
    1067                            const psVector *restrict maskVector,
    1068                            unsigned int maskVal,
    1069                            psStats *stats)
    1070 {
    1071     psHistogram *robustHistogram = NULL;
    1072     float binSize = 0.0;                        // Size of the histogram bins
    1073     float sigmaE = 0.0;
    1074     int LQBinNum = -1;                          // Bin num for lower quartile
    1075     int UQBinNum = -1;                          // Bin num for upper quartile
    1076     int i = 0;                                  // Loop index variable
    1077     int maxBinNum = 0;
    1078     int maxBinCount = 0;
    1079     float dL = 0.0;
    1080     int numBins = 0;
    1081     psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
    1082 
    1083     // NOTE: The SDRS states that the sample quartiles must be used to
    1084     // determine the initial bin sizes.  However, the sample quartiles are
    1085     // calculated based on a full sort of the data set, regardless of the
    1086     // size of the data set.  We should consult with IfA to ensure that this
    1087     // is really required.
    1088     /*
    1089         if (isnan(stats->sampleUQ) ||
    1090             isnan(stats->sampleLQ)) {
    1091             stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
    1092             p_psVectorSampleQuartiles(myVector,
    1093                                       maskVector,
    1094                                       maskVal,
    1095                                       stats);
    1096         }
    1097     */
    1098     // Compute the initial bin size of the robust histogram.
    1099     p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
    1100     binSize = stats->clippedStdev / 10.0f;
    1101 
    1102     // If stats->clippedStdev == 0.0, then all data elements have the same
    1103     // value.  Therefore, we can set the appropiate results and return.
    1104     if (fabs(binSize) <= FLT_EPSILON) {
    1105         if (stats->options & PS_STAT_ROBUST_MEAN) {
    1106             stats->robustMean = stats->clippedMean;
    1107         }
    1108         if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1109             stats->robustMedian = stats->clippedMean;
    1110         }
    1111         if  (stats->options & PS_STAT_ROBUST_MODE) {
    1112             stats->robustMode = stats->clippedMean;
    1113         }
    1114         if  (stats->options & PS_STAT_ROBUST_STDEV) {
    1115             stats->robustStdev = 0.0;
    1116         }
    1117         if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1118             stats->robustUQ = stats->clippedMean;
    1119             stats->robustLQ = stats->clippedMean;
    1120         }
    1121         psFree(tmpStats);
    1122         psFree(robustHistogram);
    1123         return;
    1124     }
    1125 
    1126     // Detemine minimum and maximum values in the data vector.
    1127     if (isnan(stats->min)) {
    1128         p_psVectorMin(myVector, maskVector, maskVal, stats);
    1129     }
    1130     if (isnan(stats->max)) {
    1131         p_psVectorMax(myVector, maskVector, maskVal, stats);
    1132     }
    1133 
    1134     // Create the histogram structure (yes, 2 is necessary, not 1).  Also,
    1135     // if we get here, we know that binSize != 0.0.
    1136     numBins = 2 + (int) ((stats->max - stats->min) / binSize);
    1137 
    1138     robustHistogram = psHistogramAlloc(stats->min,
    1139                                        stats->max,
    1140                                        numBins);
    1141     // Populate the histogram array.
    1142     psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
    1143 
    1144     // Smooth the histogram.
    1145     p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f);
    1146 
    1147     LQBinNum = -1;
    1148     UQBinNum = -1;
    1149     for (i=0;i<robustHistogram->nums->n;i++) {
    1150         if ((robustHistogram->nums->data.S32[i] <= stats->sampleLQ) &&
    1151                 (stats->sampleLQ <= robustHistogram->nums->data.S32[i])) {
    1152             LQBinNum = i;
    1153         }
    1154 
    1155         if ((robustHistogram->nums->data.S32[i] <= stats->sampleUQ) &&
    1156                 (stats->sampleUQ <= robustHistogram->nums->data.S32[i])) {
    1157             UQBinNum = i;
    1158         }
    1159     }
    1160 
    1161     // Determine the bin with the peak value in the range LQ to UQ.
    1162     maxBinNum = LQBinNum;
    1163     maxBinCount = robustHistogram->nums->data.S32[maxBinNum];
    1164     for (i=LQBinNum;i<=UQBinNum;i++) {
    1165         if (robustHistogram->nums->data.S32[i] > maxBinCount) {
    1166             maxBinNum = i;
    1167             maxBinCount = robustHistogram->nums->data.S32[i];
    1168         }
    1169     }
    1170 
    1171     dL = (stats->robustUQ - stats->robustLQ) / 8.0;
    1172 
    1173     // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    1174     // NOTE: This step is dependent on the functions in psMinimize.c being
    1175     // implemented.  Currently, they are not.
    1176 
    1177     if (stats->options & PS_STAT_ROBUST_MEAN) {
    1178         stats->robustMean = 0.0;
    1179     }
    1180 
    1181     if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1182         stats->robustMedian = 0.0;
    1183     }
    1184     if  (stats->options & PS_STAT_ROBUST_MODE) {
    1185         stats->robustMode = maxBinNum;
    1186     }
    1187     if  (stats->options & PS_STAT_ROBUST_STDEV) {
    1188         stats->robustStdev = 0.0;
    1189     }
    1190     if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1191         stats->robustUQ = 0.0;
    1192         stats->robustLQ = 0.0;
    1193     }
    1194     stats->robustNfit = 0.0;
    1195     stats->robustN50 = 0.0;
    1196     psFree(tmpStats);
    1197     psFree(robustHistogram);
    1198 }
    1199 
    12001758
    12011759/******************************************************************************
    … …  
    12291787    }
    12301788
    1231     // Ensure that the data is of type PS_TYPE_F32.  Eventually, more data
    1232     // types will be implemented.
    1233     if (in->type.type != PS_TYPE_F32) {
    1234         psAbort(__func__,
    1235                 "Only data type PS_TYPE_F32 is currently supported (0x%x).",
    1236                 in->type.type);
    1237     }
    1238 
    1239     // Ensure that the mask vector is of the proper size and type.
     1789    PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
    12401790    if (mask != NULL) {
    1241         if (in->n != mask->n) {
    1242             psAbort(__func__,
    1243                     "Vector data and vector mask are of different sizes.");
    1244         }
    1245         if (mask->type.type != PS_TYPE_U8) {
    1246             psAbort(__func__,
    1247                     "Vector mask must be type PS_TYPE_U8");
    1248         }
     1791        PS_CHECK_NULL_VECTOR(mask);
     1792        PS_CHECK_EMPTY_VECTOR(mask);
     1793        PS_CHECK_VECTOR_SIZE_EQUAL(mask, in);
     1794        PS_CHECK_VECTOR_TYPE(mask, PS_TYPE_U8);
    12491795    }
    12501796
    … …  
    12811827            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
    12821828        p_psVectorRobustStats(in, mask, maskVal, stats);
     1829
     1830        printf("HMMMM stats->robustMode is %f\n", stats->robustMode);
    12831831    }
    12841832
Note: See TracChangeset for help on using the changeset viewer.