IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1020


Ignore:
Timestamp:
Jun 14, 2004, 9:33:09 AM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r889 r1020  
    44 *         routines.
    55 *
    6  *  @author George Gusciora, MHPCC
     6 *  This file will hold the functions for allocated, freeing, and evaluating
     7 *  polynomials.  It also contains a Gaussian functions.
    78 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-06-07 00:32:53 $
     9 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-06-14 19:32:42 $
    1011 *
    1112 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    103104psPolynomial1D *psPolynomial1DAlloc(int n)
    104105{
     106    int i = 0;
    105107    psPolynomial1D *newPoly = NULL;
    106108
    … …  
    110112    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
    111113    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     114    for (i=0;i<n;i++) {
     115        newPoly->coeff[i] = 0.0;
     116        newPoly->coeffErr[i] = 0.0;
     117        newPoly->mask[i] = 0.0;
     118    }
    112119
    113120    return(newPoly);
    … …  
    117124{
    118125    int x = 0;
     126    int y = 0;
    119127    psPolynomial2D *newPoly = NULL;
    120128
    … …  
    131139        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    132140    }
    133 
     141    for (x=0;x<nX;x++) {
     142        for (y=0;y<nY;y++) {
     143            newPoly->coeff[x][y] = 0.0;
     144            newPoly->coeffErr[x][y] = 0.0;
     145            newPoly->mask[x][y] = 0;
     146        }
     147    }
    134148
    135149    return(newPoly);
    … …  
    140154    int x = 0;
    141155    int y = 0;
     156    int z = 0;
    142157    psPolynomial3D *newPoly = NULL;
    143158
    … …  
    160175        }
    161176    }
     177    for (x=0;x<nX;x++) {
     178        for (y=0;y<nY;y++) {
     179            for (z=0;z<nZ;z++) {
     180                newPoly->coeff[x][y][z] = 0.0;
     181                newPoly->coeffErr[x][y][z] = 0.0;
     182                newPoly->mask[x][y][z] = 0;
     183            }
     184        }
     185    }
    162186
    163187    return(newPoly);
    … …  
    169193    int x = 0;
    170194    int y = 0;
     195    int z = 0;
    171196    psPolynomial4D *newPoly = NULL;
    172197
    … …  
    195220        }
    196221    }
     222    for (w=0;w<nW;w++) {
     223        for (x=0;x<nX;x++) {
     224            for (y=0;y<nY;y++) {
     225                for (z=0;z<nZ;z++) {
     226                    newPoly->coeff[w][x][y][z] = 0.0;
     227                    newPoly->coeffErr[w][x][y][z] = 0.0;
     228                    newPoly->mask[w][x][y][z] = 0;
     229                }
     230            }
     231        }
     232    }
    197233
    198234    return(newPoly);
    … …  
    394430psDPolynomial1D *psDPolynomial1DAlloc(int n)
    395431{
     432    int i = 0;
    396433    psDPolynomial1D *newPoly = NULL;
    397434
    … …  
    401438    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
    402439    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     440    for (i=0;i<n;i++) {
     441        newPoly->coeff[i] = 0.0;
     442        newPoly->coeffErr[i] = 0.0;
     443        newPoly->mask[i] = 0.0;
     444    }
    403445
    404446    return(newPoly);
    … …  
    408450{
    409451    int x = 0;
     452    int y = 0;
    410453    psDPolynomial2D *newPoly = NULL;
    411454
    … …  
    422465        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    423466    }
     467    for (x=0;x<nX;x++) {
     468        for (y=0;y<nY;y++) {
     469            newPoly->coeff[x][y] = 0.0;
     470            newPoly->coeffErr[x][y] = 0.0;
     471            newPoly->mask[x][y] = 0;
     472        }
     473    }
    424474
    425475    return(newPoly);
    … …  
    430480    int x = 0;
    431481    int y = 0;
     482    int z = 0;
    432483    psDPolynomial3D *newPoly = NULL;
    433484
    … …  
    450501        }
    451502    }
     503    for (x=0;x<nX;x++) {
     504        for (y=0;y<nY;y++) {
     505            for (z=0;z<nZ;z++) {
     506                newPoly->coeff[x][y][z] = 0.0;
     507                newPoly->coeffErr[x][y][z] = 0.0;
     508                newPoly->mask[x][y][z] = 0;
     509            }
     510        }
     511    }
    452512
    453513    return(newPoly);
    … …  
    459519    int x = 0;
    460520    int y = 0;
     521    int z = 0;
    461522    psDPolynomial4D *newPoly = NULL;
    462523
    … …  
    485546        }
    486547    }
     548    for (w=0;w<nW;w++) {
     549        for (x=0;x<nX;x++) {
     550            for (y=0;y<nY;y++) {
     551                for (z=0;z<nZ;z++) {
     552                    newPoly->coeff[w][x][y][z] = 0.0;
     553                    newPoly->coeffErr[w][x][y][z] = 0.0;
     554                    newPoly->mask[w][x][y][z] = 0;
     555                }
     556            }
     557        }
     558    }
    487559
    488560    return(newPoly);

  • trunk/psLib/src/dataManip/psFunctions.h

    r974 r1020  
     1/** @file psFunctions.h
     2 *  \brief Standard Mathematical Functions.
     3 *  \ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author George Gusciora, MHPCC
     13 *
     14 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2004-06-14 19:32:42 $
     16 *
     17 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     18 */
     19
    120#if !defined(PS_FUNCTIONS_H)
    221#define PS_FUNCTIONS_H
    322
    4 /** \file psFunctions.h
    5  *  \brief Standard Mathematical Functions.
    6  *  \ingroup Stats
    7  */
    823
    924/** \addtogroup Stats
    … …  
    2338
    2439
    25 psVector *psGaussianDev(float mean,
    26                         float sigma,
    27                         int Npts);
    28 
    29 /*****************************************************************************/
     40/** Produce a vector of random numbers from a Gaussian distribution with
     41    the specified mean and sigma */
     42psVector *psGaussianDev(float mean,    ///< The mean of the Gaussian
     43                        float sigma,   ///< The sigma of the Gaussian
     44                        int Npts);     ///< The size of the vector
     45
     46
     47
     48
    3049
    3150/** One-dimensional polynomial */
    3251typedef struct
    3352{
    34     int n;    ///< Number of terms
    35     float *coeff;   ///< Coefficients
    36     float *coeffErr;   ///< Error in coefficients
    37     char *mask;    ///< Coefficient mask
     53    int n;           ///< Number of terms
     54    float *coeff;    ///< Coefficients
     55    float *coeffErr; ///< Error in coefficients
     56    char *mask;      ///< Coefficient mask
    3857}
    3958psPolynomial1D;

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

    r978 r1020  
    1515
    1616#include "float.h"
    17 #define ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
     17#define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
    1818// will use robust statistical methods.
    19 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar
    20 // smoothing.
    21 #define PI 3.141592653
     19#define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
    2220#define CLIPPED_NUM_ITER_LB 1
    2321#define CLIPPED_NUM_ITER_UB 10
    … …  
    2624#define true 1
    2725#define false 0
    28 #define MYMAXFLOAT 1e99
     26#define MYMAXFLOAT HUGE
    2927
    3028#ifndef DOXYGEN
    … …  
    4846    newStruct->sampleUQ = NAN;
    4947    newStruct->sampleLQ = NAN;
    50     newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD;
     48    newStruct->sampleLimit = DEFAULT_ROBUST_SIZE_THRESHOLD;
    5149    newStruct->robustMean = NAN;
    5250    newStruct->robustMedian = NAN;
    … …  
    7876
    7977/******************************************************************************
    80     HISTOGRAM FUNCTIONS
     78psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure
     79with the specifed upper and lower limits, and the specifed number of bins.
     80This routine will also set the bounds for each of the bins.
     81 
     82Input:
     83    lower
     84    upper
     85    n
     86Returns:
     87    The histogram structure
    8188 *****************************************************************************/
    8289psHistogram *psHistogramAlloc(float lower,
    … …  
    8491                              int n)
    8592{
    86     int i = 0;
    87     psHistogram *newHist = NULL;
    88     float binSize = 0.0;
    89 
     93    int i = 0;                          // Loop index variable
     94    psHistogram *newHist = NULL;        // The new histogram structure
     95    float binSize = 0.0;                // The histogram bin size
     96
     97    // GUS: Verify that this is the correct action.
    9098    if (n == 0) {
    91         psAbort(__func__, "psHistogram requested with 0 bins");
     99        return(NULL);
    92100    }
    93101
    94102    // GUS: Verify that this is the correct action.
    95103    if (lower > upper) {
    96         //        psAbort(__func__, "psHistogram requested with lower>upper");
    97104        return(NULL);
    98105    }
    99106
     107    // Allocate memory for the new histogram structure.  If there are N
     108    // bins, then there are N+1 bounds to those bins.
    100109    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    101110    newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32);
    102111    newHist->bounds->n = newHist->bounds->nalloc;
    103112
     113    // Calculate the bounds for each bin.
    104114    binSize = (upper - lower) / (float) n;
    105115    for (i=0;i<n+1;i++) {
    106116        newHist->bounds->data.F32[i] = lower + (binSize * (float) i);
    107117    }
     118
     119    // Allocate the bins, and initialize them to zero.
    108120    newHist->nums =  psVectorAlloc(n, PS_TYPE_S32);
    109121    newHist->nums->n = newHist->nums->nalloc;
    … …  
    111123        newHist->nums->data.S32[i] = 0;
    112124    }
     125
     126    // Initialize the other members.
    113127    newHist->minNum = 0;
    114128    newHist->maxNum = 0;
    … …  
    118132}
    119133
    120 // When this is called, the number of data elements in bounds
    121 // is n.  Therefore, the number of bins is n-1.
     134/******************************************************************************
     135psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure
     136with the specifed bounds.  The number of elements in the bounds vector is n.
     137Therefore, the number of bins is n-1.
     138Input:
     139    bounds
     140Returns:
     141    The histogram structure
     142 *****************************************************************************/
    122143psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
    123144{
    124     psHistogram *newHist = NULL;
    125     int i;
     145    psHistogram *newHist = NULL;        // The new histogram structure
     146    int i;                              // Loop index variable
    126147
    127148    // GUS: Verify that this is the correct action.
    … …  
    131152    }
    132153
     154    // GUS: Verify that this is the correct action.
     155    if (bounds->n <= 1) {
     156        //        psAbort(__func__, "psHistogram requested with NULL bounds");
     157        return(NULL);
     158    }
     159
     160    // Allocate memory for the new histogram structure.
    133161    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    134162    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
    … …  
    137165        newHist->bounds->data.F32[i] = bounds->data.F32[i];
    138166    }
     167
     168    // Allocate the bins, and initialize them to zero.  If there are N bounds,
     169    // then there are N-1 bins.
    139170    newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32);
    140171    newHist->nums->n = newHist->nums->nalloc;
    … …  
    143174    }
    144175
     176    // Initialize the other members.
    145177    newHist->minNum = 0;
    146178    newHist->maxNum = 0;
    … …  
    159191
    160192/*****************************************************************************
    161     NOTE: Can we assume fixed size bins in the myHist data structure?
     193psHistogramVector(out, in, mask, maskVal): this procedure takes as input a
     194preallocated and initialized histogram structure.  It fills the bins in that
     195histogram structure in accordance with the input data "in" and the, possibly
     196NULL, mask vector.
    162197 
    163     NOTE: This procedure assumes that the following has been already set
    164         myHist->lower
    165         myHist->upper
    166  myHist->minVal
    167  myHist->maxVal
    168     This procedure sets
    169         myHist->nums
    170         myHist->minNum
    171         myHist->maxNum
     198Inputs:
     199    out
     200    in
     201    mask
     202    maskVal
     203Returns:
     204    The histogram structure "out".
    172205 *****************************************************************************/
    173206psHistogram *psHistogramVector(psHistogram *out,
    174                                psVector *in,
    175                                psVector *mask,
     207                               const psVector *restrict in,
     208                               const psVector *restrict mask,
    176209                               int maskVal)
    177210{
    178     int i = 0;
    179     int j = 0;
    180     float binSize = 0.0;
    181     int binNum = 0;
    182     int numBins = 0;
     211    int i = 0;                                  // Loop index variable
     212    int j = 0;                                  // Loop index variable
     213    float binSize = 0.0;                        // Histogram bin size
     214    int binNum = 0;                             // A temporary bin number
     215    int numBins = 0;                            // The total number of bins
    183216
    184217    // GUS: Verify that this is the correct action.
    185218    if (out == NULL) {
    186         //        psAbort(__func__, "psHistogram requested with NULL struct");
    187219        return(NULL);
    188220    }
    … …  
    190222    // GUS: Verify that this is the correct action.
    191223    if (in == NULL) {
    192         //        psAbort(__func__, "psHistogram requested with NULL struct");
    193224        return(out);
    194225    }
    195226
     227    if (in->type.type != PS_TYPE_F32) {
     228        psAbort(__func__,
     229                "Only data type PS_TYPE_F32 is currently supported (0x%x).",
     230                in->type.type);
     231    }
     232
     233    if (mask != NULL) {
     234        if (in->n != mask->n) {
     235            psAbort(__func__,
     236                    "Vector data and vector mask are of different sizes.");
     237        }
     238        if (mask->type.type != PS_TYPE_U8) {
     239            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
     240        }
     241    }
     242    // GUS: determine the correct action for a variety of other cases:
     243    // in vector has 0 elements, and histogram structure has zero bins.
     244
    196245    numBins = out->nums->n;
    197246
    198     if (mask != NULL) {
    199         for (i=0;i<in->n;i++) {
    200             // Check if this pixel is masked, and if so, skip it.
    201             if (!(mask->data.U8[i] & maskVal)) {
    202                 // Check if this pixel is below the minimum value, and if so
    203                 // count it, then skip it.
    204                 if (in->data.F32[i] < out->bounds->data.F32[0]) {
    205                     out->minNum++;
    206 
    207                     // Check if this pixel is above the maximum value, and if so
    208                     // count it, then skip it.
    209                 } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    210                     out->maxNum++;
    211                 } else {
    212                     // If this is a uniform histogram, determining the correct
    213                     // number is trivial.
    214                     if (out->uniform == true) {
    215                         binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    216 
    217                         binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    218                                         binSize);
    219                         (out->nums->data.S32[binNum])++;
    220                         // If this is a non-uniform histogram, determining the correct
    221                         // bin number requires a bit more work.
    222                     } else {
    223                         // GUS: This is slow.  Put a smarter algorithm here to
    224                         // find the correct bin number (bin search, probably)
    225                         for (j=0;j<(out->bounds->n)-1;j++) {
    226                             if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
    227                                     (in->data.F32[i] <= out->bounds->data.S32[j+1])) {
    228                                 (out->nums->data.S32[j])++;
    229                             }
    230                         }
    231                     }
    232                 }
    233             }
    234         }
    235     } else {
    236         for (i=0;i<in->n;i++) {
     247    for (i=0;i<in->n;i++) {
     248        // Check if this pixel is masked, and if so, skip it.
     249        if ((mask == NULL) ||
     250                ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
     251            // Check if this pixel is below the minimum value, and if so
     252            // count it, then skip it.
    237253            if (in->data.F32[i] < out->bounds->data.F32[0]) {
    238                 // Check if this pixel is below the minimum value, and if so
    239                 // count it, then skip it.
    240254                out->minNum++;
    241             } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    242255                // Check if this pixel is above the maximum value, and if so
    243256                // count it, then skip it.
     257            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    244258                out->maxNum++;
    245259            } else {
    … …  
    251265                    binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    252266                                    binSize);
    253                     (out->nums->data.S32[binNum])++;
    254 
    255267                    // If this is a non-uniform histogram, determining the correct
    256268                    // bin number requires a bit more work.
    … …  
    267279            }
    268280        }
    269 
    270281    }
    271282    return(out);
    272283}
    273284
    274 void p_printVector(psVector *myVector,
    275                    psVector *maskVector,
    276                    unsigned int maskVal,
    277                    psStats *stats)
    278 {
    279     int i = 0;
     285/*****************************************************************************
     286p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
     287function that simply prints a vector to STDOUT.  Used primarily for
     288debugging.
     289 *****************************************************************************/
     290
     291void p_psVectorPrint(psVector *myVector,
     292                     psVector *maskVector,
     293                     unsigned int maskVal,
     294                     psStats *stats)
     295{
     296    int i = 0;                                  // Loop index variable.
    280297
    281298    for (i=0;i<myVector->n;i++) {
    … …  
    291308 ******************************************************************************
    292309 ******************************************************************************
    293     MISC STATISTICAL FUNCTIONS
     310    MISC PRIVATE STATISTICAL FUNCTIONS
     311 
     312    NOTE: it is assumed that any call to these statistical functions will
     313    have been preceded by a call to the psVectorStats() function.  Various
     314    sanity tests will only be performed in psVectorStats().
     315        Is the mask vector the same length as the data vector?
     316        Is the mask vector of type PS_TYPE_U8?
     317        Is the stats data structure NULL?
     318        Is the in data structure NULL?
     319        Is the in data structure of type PS_TYPE_F32?
    294320 ******************************************************************************
    295321 ******************************************************************************
    … …  
    298324
    299325/******************************************************************************
    300     ASSUMPTION: the mean is always calculated exactly.  Robust means are never
    301     calculated in this routine.
     326p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
     327mean of the input vector.
     328Inputs
     329Returns
     330    NULL
     331ASSUMPTION: the mean is always calculated exactly.  Robust means are never
     332calculated in this routine.
    302333 *****************************************************************************/
    303334void p_psVectorSampleMean(const psVector *restrict myVector,
    … …  
    306337                          psStats *stats)
    307338{
    308     int i = 0;
    309     float mean = 0.0;
    310     int count = 0;
    311     float rangeMin = 0.0;
    312     float rangeMax = 0.0;
    313 
     339    int i = 0;                                  // Loop index variable
     340    float mean = 0.0;                           // The mean
     341    int count = 0;                              // # of points in this mean?
     342    float rangeMin = 0.0;                       // Exclude data below this
     343    float rangeMax = 0.0;                       // Exclude date above this
     344
     345    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
     346    // loop.
    314347    if (stats->options & PS_STAT_USE_RANGE) {
    315348        rangeMin = stats->min;
    … …  
    317350        if (maskVector != NULL) {
    318351            for (i=0;i<myVector->n;i++) {
     352                // Check if the data is with the specified range
    319353                if (!(maskVal & maskVector->data.U8[i]) &&
    320354                        (rangeMin <= myVector->data.F32[i]) &&
    … …  
    355389}
    356390
     391/******************************************************************************
     392p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     393max of the input vector.
     394Inputs
     395    myVector
     396    maskVector
     397    maskVal
     398    stats
     399Returns
     400    NULL
     401 *****************************************************************************/
    357402void p_psVectorMax(const psVector *restrict myVector,
    358403                   const psVector *restrict maskVector,
    … …  
    360405                   psStats *stats)
    361406{
    362     int i = 0;
    363     float max = -MYMAXFLOAT;
    364     float rangeMin = 0.0;
    365     float rangeMax = 0.0;
    366 
     407    int i = 0;                                  // Loop index variable
     408    float max = -MYMAXFLOAT;                    // The calculated maximum
     409    float rangeMin = 0.0;                       // Exclude data below this
     410    float rangeMax = 0.0;                       // Exclude date above this
     411
     412    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    367413    if (stats->options & PS_STAT_USE_RANGE) {
    368414        rangeMin = stats->min;
    … …  
    408454}
    409455
     456/******************************************************************************
     457p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
     458minimum of the input vector.
     459Inputs
     460    myVector
     461    maskVector
     462    maskVal
     463    stats
     464Returns
     465    NULL
     466 *****************************************************************************/
    410467void p_psVectorMin(const psVector *restrict myVector,
    411468                   const psVector *restrict maskVector,
    … …  
    413470                   psStats *stats)
    414471{
    415     int i = 0;
    416     float min = MYMAXFLOAT;
    417     float rangeMin = 0.0;
    418     float rangeMax = 0.0;
     472    int i = 0;                                  // Loop index variable
     473    float min = MYMAXFLOAT;                     // The calculated maximum
     474    float rangeMin = 0.0;                       // Exclude data below this
     475    float rangeMax = 0.0;                       // Exclude date above this
    419476
    420477    if (stats->options & PS_STAT_USE_RANGE) {
    … …  
    462519
    463520/******************************************************************************
    464     This routine calculates the number of non-masked pixels in the vector.
    465     that fall within the min/max range, if given.
     521p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     522number of non-masked pixels in the vector.  that fall within the min/max
     523range, if given.
     524Inputs
     525    myVector
     526    maskVector
     527    maskVal
     528    stats
     529Returns
     530    NULL
    466531 *****************************************************************************/
    467532int p_psVectorNValues(const psVector *restrict myVector,
    … …  
    470535                      psStats *stats)
    471536{
    472     int i = 0;
    473     int numData = 0;
    474     float rangeMin = 0.0;
    475     float rangeMax = 0.0;
     537    int i = 0;                                  // Loop index variable
     538    int numData = 0;                            // The number of data points
     539    float rangeMin = 0.0;                       // Exclude data below this
     540    float rangeMax = 0.0;                       // Exclude date above this
    476541
    477542    if (stats->options & PS_STAT_USE_RANGE) {
    … …  
    512577
    513578
     579/******************************************************************************
     580p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
     581median of the input vector.
     582Inputs
     583    myVector
     584    maskVector
     585    maskVal
     586    stats
     587Returns
     588    NULL
     589 *****************************************************************************/
    514590void p_psVectorSampleMedian(const psVector *restrict myVector,
    515591                            const psVector *restrict maskVector,
    … …  
    517593                            psStats *stats)
    518594{
    519     psVector *unsortedVector = NULL;
    520     psVector *sortedVector = NULL;
    521     int count = 0;
    522     int i = 0;
    523     int nValues = 0;
    524     float rangeMin = 0.0;
    525     float rangeMax = 0.0;
    526     psStats *stats2 = NULL;
     595    psVector *unsortedVector = NULL;            // Temporary vector
     596    psVector *sortedVector = NULL;              // Temporary vector
     597    int i = 0;                                  // Loop index variable
     598    int count = 0;                              // # of points in this mean?
     599    int nValues = 0;                            // # of points in vector
     600    float rangeMin = 0.0;                       // Exclude data below this
     601    float rangeMax = 0.0;                       // Exclude date above this
     602    psStats *stats2 = NULL;                     // Temporary stats structure
    527603
    528604
    … …  
    531607
    532608    if (myVector->n > stats->sampleLimit) {
     609        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     610
    533611        // Calculate the robust quartiles.
    534612        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    … …  
    548626
    549627    // Determine how many data points fit inside this min/max range
    550     // and are not maxed, IF the maskVector is not NULL>
     628    // and are not masked, IF the maskVector is not NULL>
    551629    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    552630
    … …  
    620698                                  float sigma)
    621699{
    622     int i = 0;
    623     int j = 0;
    624     float tmpf = 0.0;
    625     float gaussianCoefs[1 + (2 * GAUSS_WIDTH)];
    626     // The coefficients used in the histogram
    627     // smoothing calculation.
     700    int i = 0;                                  // Loop index variable
     701    int j = 0;                                  // Loop index variable
     702    float tmpf = 0.0;                           // Temporary variable
     703    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
    628704
    629705    for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
    … …  
    631707            // If sigma does not equal zero, then we use Gaussian smoothing.
    632708            #ifdef  DARWIN
    633             tmpf = (float) sqrt(2.0f * PI * sigma * sigma);
     709            tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma);
    634710            #else
    635711
    636             tmpf = sqrtf(2.0f * PI * sigma * sigma);
     712            tmpf = sqrtf(2.0f * M_PI * sigma * sigma);
    637713            #endif
    638714
    639715            gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) *
    640716                                             ((float) (i-GAUSS_WIDTH))) /
    641                                             (2.0f * sigma * sigma)) /
    642                                tmpf;
     717                                            (2.0f * sigma * sigma)) / tmpf;
    643718        } else {
    644719            /* If sigma equals zero (all pixels have the same value)
    645              * the above code will divide by zero.  Therefore, we instead
    646              * use boxcar smoothing.
     720             * the above code will divide by zero.  Therefore, we don't need
     721             * to smooth the data.
    647722             */
    648             gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH));
     723            return;
    649724        }
    650725    }
    … …  
    662737
    663738/******************************************************************************
    664     p_psVectorSampleQuartiles()
    665  This procedure calculates the upper and/or lower quartiles of the
    666  data set.
    667  
     739p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
     740the upper and/or lower quartiles of the input vector.
     741Inputs
     742    myVector
     743    maskVector
     744    maskVal
     745    stats
     746Returns
     747    NULL
    668748 *****************************************************************************/
    669749void p_psVectorSampleQuartiles(const psVector *restrict myVector,
    … …  
    672752                               psStats *stats)
    673753{
    674     psVector *unsortedVector = NULL;
    675     psVector *sortedVector = NULL;
    676     int count = 0;
    677     int ind = 0;
    678     int i = 0;
    679     int nValues = 0;
    680     float rangeMin = 0.0;
    681     float rangeMax = 0.0;
     754    psVector *unsortedVector = NULL;            // Temporary vector
     755    psVector *sortedVector = NULL;              // Temporary vector
     756    int i = 0;                                  // Loop index variable
     757    int count = 0;                              // # of points in this mean?
     758    int nValues = 0;                            // # data points
     759    float rangeMin = 0.0;                       // Exclude data below this
     760    float rangeMax = 0.0;                       // Exclude date above this
    682761
    683762    // Determine if the number of data points exceed a threshold which will
    684763    // cause to generate robust stats, as opposed to exact stats.
    685     /* GUS: When IfA provides the algorithm, insert is here.
    686         if (myVector->n > stats->sampleLimit) {
    687             psStats *stats2 = NULL;
    688             // Calculate the robust quartiles.
    689             stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
    690             p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    691      
    692             // Store the robust quartiles into the sample quartile members.
    693             stats->sampleUQ = stats2->robustUQ;
    694             stats->sampleLQ = stats2->robustLQ;
    695      
    696             // Free temporary data buffers.
    697             psStatsFree(stats2);
    698      
    699             // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    700             stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    701      
    702             return;
    703         }
    704     */
     764    if (myVector->n > stats->sampleLimit) {
     765        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     766        psStats *stats2 = NULL;
     767        // Calculate the robust quartiles.
     768        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
     769        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     770
     771        // Store the robust quartiles into the sample quartile members.
     772        stats->sampleUQ = stats2->robustUQ;
     773        stats->sampleLQ = stats2->robustLQ;
     774
     775        // Free temporary data buffers.
     776        psStatsFree(stats2);
     777
     778        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     779        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     780
     781        return;
     782    }
    705783
    706784    // Determine how many data points fit inside this min/max range
    … …  
    757835
    758836    // Calculate the quartile points exactly.
    759     ind = 3 * (nValues / 4);
    760     stats->sampleUQ = sortedVector->data.F32[ind];
    761     ind = (nValues / 4);
    762     stats->sampleLQ = sortedVector->data.F32[ind];
     837    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     838    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    763839
    764840    // Free the temporary data structures.
    … …  
    770846
    771847/******************************************************************************
    772     p_psVectorRobustStats(): this procedure calculates a variety of robust
    773     stat measures:
     848p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
     849calculates a variety of robust stat measures:
    774850        PS_STAT_ROBUST_MEAN
    775851        PS_STAT_ROBUST_MEDIAN
    … …  
    784860    measures were calcualted in separate functiosn, then much of the initial
    785861    processing would be duplicated.
     862Inputs
     863    myVector
     864    maskVector
     865    maskVal
     866    stats
     867Returns
     868    NULL
    786869 *****************************************************************************/
    787870void p_psVectorRobustStats(const psVector *restrict myVector,
    … …  
    791874{
    792875    psHistogram *robustHistogram = NULL;
    793     float binSize = 0.0;
     876    float binSize = 0.0;                        // Size of the histogram bins
    794877    float sigmaE = 0.0;
    795     int LQBinNum = -1;
    796     int UQBinNum = -1;
    797     int i = 0;    // Loop index variable.
     878    int LQBinNum = -1;                          // Bin num for lower quartile
     879    int UQBinNum = -1;                          // Bin num for upper quartile
     880    int i = 0;                                  // Loop index variable
    798881    int maxBinNum = 0;
    799882    int maxBinCount = 0;
    800883    float dL = 0.0;
     884    int numBins = 0;
    801885
    802886    // NOTE: The SDRS states that the sample quartiles must be used to
    … …  
    827911    }
    828912
    829     // Create the histogram structure.
     913    // Create the histogram structure (yes, 2 is necessary, not 1).
     914    numBins = 2 + (int) ((stats->max - stats->min) / binSize);
     915
    830916    robustHistogram = psHistogramAlloc(stats->min,
    831917                                       stats->max,
    832                                        1 + (int) ((stats->max - stats->min) / binSize));
    833 
    834 
     918                                       numBins);
    835919    // Populate the histogram array.
    836     // GUS: fix this
    837     //    robustHistogram = psHistogramVector(robustHistogram, myVector);
    838     //
     920    psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
    839921
    840922    // Smooth the histogram.
    … …  
    868950
    869951    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    870     // What algorithm should I use to do this?
    871 
    872     // The resulting fit parameters are the robust mean, mean_r, and sigma
    873     // What is the mean_r?
     952    // GUS: This step is dependent on the functions in psMinimize.c being
     953    // implemented.  Currently, they are not.
    874954
    875955    if (stats->options & PS_STAT_ROBUST_MEAN) {
    … …  
    895975}
    896976
    897 /*****************************************************************************
    898     NOTE: This function assumes that p_psVectorMean() has already been called
    899     and the correct value is stored in stats->sampleMean.
     977/******************************************************************************
     978p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
     979stdev of the input vector.
     980Inputs
     981    myVector
     982    maskVector
     983    maskVal
     984    stats
     985Returns
     986    NULL
    900987 
    901     NOTE: the mean is always calculated exactly.  Robust means are never
    902     calculated in this routine.
     988NOTE: the mean is always calculated exactly.  Robust means are never
     989calculated in this routine.
    903990 *****************************************************************************/
    904991void p_psVectorSampleStdev(const psVector *restrict myVector,
    … …  
    907994                           psStats *stats)
    908995{
    909     int i            = 0;
    910     int countInt     = 0;
    911     float countFloat = 0.0;
    912     float       mean = 0.0;
    913     float       diff = 0.0;
    914     float sumSquares = 0.0;
    915     float sumDiffs = 0.0;
    916     float rangeMin = 0.0;
    917     float rangeMax = 0.0;
    918 
     996    int i = 0;                                  // Loop index variable
     997    int countInt = 0;                           // # of data points being used
     998    float countFloat = 0.0;                     // # of data points being used
     999    float mean = 0.0;                           // The mean
     1000    float diff = 0.0;                           // Used in calculating stdev
     1001    float sumSquares = 0.0;                     // temporary variable
     1002    float sumDiffs = 0.0;                       // temporary variable
     1003    float rangeMin = 0.0;                       // Exclude data below this
     1004    float rangeMax = 0.0;                       // Exclude date above this
     1005
     1006    // This procedure requires the mean.  If it has not been already
     1007    // calculated, then call p_psVectorSampleMean()
    9191008    if (0 != isnan(stats->sampleMean)) {
    9201009        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    … …  
    9791068}
    9801069
     1070/******************************************************************************
     1071p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
     1072clipped stats (mean or stdev) of the input vector.
     1073 
     1074Inputs
     1075    myVector
     1076    maskVector
     1077    maskVal
     1078    stats
     1079Returns
     1080    NULL
     1081 *****************************************************************************/
    9811082void p_psVectorClippedStats(const psVector *restrict myVector,
    9821083                            const psVector *restrict maskVector,
    … …  
    9841085                            psStats *stats)
    9851086{
    986     int i = 0;
    987     int j = 0;
    988     float clippedMean = 0.0;
    989     float clippedStdev = 0.0;
    990     psVector *tmpMask = NULL;
    991 
     1087    int i = 0;                                  // Loop index variable
     1088    int j = 0;                                  // Loop index variable
     1089    float clippedMean = 0.0;                    // self-explanatory
     1090    float clippedStdev = 0.0;                   // self-explanatory
     1091    float oldStanMean = 0.0;                    // Temporary variable
     1092    float oldStanStdev = 0.0;                   // Temporary variable
     1093    psVector *tmpMask = NULL;                   // Temporary vector
     1094
     1095    // Endure that stats->clipIter is within the proper range.
    9921096    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
    9931097            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
    … …  
    9961100    }
    9971101
     1102    // Endure that stats->clipSigma is within the proper range.
    9981103    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
    9991104            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
    … …  
    10021107    }
    10031108
    1004 
     1109    // We allocate a temporary mask vector since during the iterative
     1110    // steps that follow, we will be masking off additional data points.
     1111    // However, we do no want to modify the original mask vector.
    10051112    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    10061113    tmpMask->n = myVector->n;
    10071114
     1115    // If we were called with a mask vector, then initialize the temporary
     1116    // mask vector with those values.
    10081117    if (maskVector != NULL) {
    10091118        for (i=0;i<tmpMask->n;i++) {
    … …  
    10261135    clippedStdev = stats->sampleStdev;
    10271136
     1137    // Must save the old sampleMean and sampleStdev since the following code
     1138    // block overwrites them.
     1139    oldStanMean = stats->sampleMean;
     1140    oldStanStdev = stats->sampleStdev;
     1141
    10281142    // 5. Repeat N times:
    1029 
    10301143    for (i=0;i<stats->clipIter;i++) {
    1031         //        printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i,
    1032         //                clippedMean, clippedStdev);
    10331144        for (j=0;j<myVector->n;j++) {
    10341145            // a) Exclude all values x_i for which |x_i - x| > K * stdev
    … …  
    10381149            }
    10391150            // b) compute new mean and stdev
    1040             // GUS: I should probably create a new struct here since the
    1041             // following calls will overwrite any old values in sampleMean.
    10421151            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
    10431152            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
    … …  
    10501159        }
    10511160    }
     1161    stats->sampleMean = oldStanMean;
     1162    stats->sampleStdev= oldStanStdev;
    10521163
    10531164    // 7. The last calcuated value of x is the cliped mean.
    … …  
    10661177
    10671178/******************************************************************************
     1179psVectorStats(myVector, maskVector, maskVal, stats): this is the public API
     1180function which calls the above private stats functions based on what bits
     1181were set in stats->options.
    10681182 
    1069     NOTE: The current strategy is to implement everything assuming that all
    1070     input data is of type PS_TYPE_F32.  Once the basic code is in place,
    1071     we will macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
     1183Inputs
     1184    myVector
     1185    maskVector
     1186    maskVal
     1187    stats
     1188Returns
     1189    The stats structure.
    10721190 
     1191NOTE: The current strategy is to implement everything assuming that all
     1192input data is of type PS_TYPE_F32.  Once the basic code is in place, we will
     1193macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
    10731194 *****************************************************************************/
    10741195psStats *psVectorStats(psStats *stats,
    … …  
    10771198                       unsigned int maskVal)
    10781199{
     1200    // GUS: Verify that this is the correct action.
    10791201    if (in == NULL) {
    10801202        return(stats);
    … …  
    10841206    }
    10851207
     1208    // Ensure that the data is of type PS_TYPE_F32.  Eventually, more data
     1209    // types will be implemented.
    10861210    if (in->type.type != PS_TYPE_F32) {
    10871211        psAbort(__func__,
    … …  
    10901214    }
    10911215
     1216    // Ensure that the mask vector is of the proper size and type.
    10921217    if (mask != NULL) {
    10931218        if (in->n != mask->n) {
    … …  
    11241249    }
    11251250
     1251    // Since the various robust stats quantities share much computation, they
     1252    // are grouped together in a single private function:
     1253    // p_psVectorRobustStats()
    11261254    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
    11271255            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
    … …  
    11471275    }
    11481276
    1149     //    OLD CODE: Should we check for an unknown option?
    1150     //    default:
    1151     //        psAbort(__func__, "Unknown options 0x%x.\n", stats->options);
    1152 
    11531277    return(stats);
    11541278}

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

    r974 r1020  
     1/** @file  psStats.h
     2 *  \brief basic statistical operations
     3 *  @ingroup Stats
     4 *
     5 *  This file will hold the definition of the histogram and stats data
     6 *  structures.  It also contains prototypes for procedures which operate
     7 *  on those data structures.
     8 *
     9 *  @author George Gusciora, MHPCC
     10 *
     11 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-06-14 19:33:09 $
     13 *
     14 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     15 */
    116#if !defined(PS_STATS_H)
    217#define PS_STATS_H
    3 /******************************************************************************
    4     This file will histogram/statistical functions and data structures.
    5  *****************************************************************************/
    6 /** \file psStats.h
    7  *  \brief basic statistical operations
    8  *  \ingroup Stats
    9  */
     18
    1019#include "psVector.h"
    1120
    … …  
    3645
    3746
    38 /** This is the generic statistics structure */
     47/** This is the generic statistics structure.  It contails the data members
     48    for the various statistic values.  It contains the options member to
     49    specifiy which statistics must be calculated. */
    3950typedef struct
    4051{
    … …  
    7485             );
    7586
    76 /** Constructor */
     87/** A constructor for the stats structure.*/
    7788psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure
    7889
    79 /** Destructor */
     90/** A destructor for the stats structure.*/
    8091void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy
    8192
    … …  
    8394    Histogram functions and data structures.
    8495 *****************************************************************************/
     96/** The basic histogram structure which contains bounds and bins. */
    8597typedef struct
    8698{
    … …  
    111123/** Calculate a histogram \ingroup MathGroup **/
    112124psHistogram *psHistogramVector (psHistogram *out,   ///< Histogram data
    113                                 psVector *in,       ///< Vector to analyse
    114                                 psVector *mask,     ///< Mask dat for input vector
     125                                const psVector *restrict in,       ///< Vector to analyse
     126                                const psVector *restrict mask,     ///< Mask dat for input vector
    115127                                int maskVal);       ///< Mask value
    116128

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

    r889 r1020  
    44 *         routines.
    55 *
    6  *  @author George Gusciora, MHPCC
     6 *  This file will hold the functions for allocated, freeing, and evaluating
     7 *  polynomials.  It also contains a Gaussian functions.
    78 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-06-07 00:32:53 $
     9 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-06-14 19:32:42 $
    1011 *
    1112 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    103104psPolynomial1D *psPolynomial1DAlloc(int n)
    104105{
     106    int i = 0;
    105107    psPolynomial1D *newPoly = NULL;
    106108
    … …  
    110112    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
    111113    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     114    for (i=0;i<n;i++) {
     115        newPoly->coeff[i] = 0.0;
     116        newPoly->coeffErr[i] = 0.0;
     117        newPoly->mask[i] = 0.0;
     118    }
    112119
    113120    return(newPoly);
    … …  
    117124{
    118125    int x = 0;
     126    int y = 0;
    119127    psPolynomial2D *newPoly = NULL;
    120128
    … …  
    131139        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    132140    }
    133 
     141    for (x=0;x<nX;x++) {
     142        for (y=0;y<nY;y++) {
     143            newPoly->coeff[x][y] = 0.0;
     144            newPoly->coeffErr[x][y] = 0.0;
     145            newPoly->mask[x][y] = 0;
     146        }
     147    }
    134148
    135149    return(newPoly);
    … …  
    140154    int x = 0;
    141155    int y = 0;
     156    int z = 0;
    142157    psPolynomial3D *newPoly = NULL;
    143158
    … …  
    160175        }
    161176    }
     177    for (x=0;x<nX;x++) {
     178        for (y=0;y<nY;y++) {
     179            for (z=0;z<nZ;z++) {
     180                newPoly->coeff[x][y][z] = 0.0;
     181                newPoly->coeffErr[x][y][z] = 0.0;
     182                newPoly->mask[x][y][z] = 0;
     183            }
     184        }
     185    }
    162186
    163187    return(newPoly);
    … …  
    169193    int x = 0;
    170194    int y = 0;
     195    int z = 0;
    171196    psPolynomial4D *newPoly = NULL;
    172197
    … …  
    195220        }
    196221    }
     222    for (w=0;w<nW;w++) {
     223        for (x=0;x<nX;x++) {
     224            for (y=0;y<nY;y++) {
     225                for (z=0;z<nZ;z++) {
     226                    newPoly->coeff[w][x][y][z] = 0.0;
     227                    newPoly->coeffErr[w][x][y][z] = 0.0;
     228                    newPoly->mask[w][x][y][z] = 0;
     229                }
     230            }
     231        }
     232    }
    197233
    198234    return(newPoly);
    … …  
    394430psDPolynomial1D *psDPolynomial1DAlloc(int n)
    395431{
     432    int i = 0;
    396433    psDPolynomial1D *newPoly = NULL;
    397434
    … …  
    401438    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
    402439    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     440    for (i=0;i<n;i++) {
     441        newPoly->coeff[i] = 0.0;
     442        newPoly->coeffErr[i] = 0.0;
     443        newPoly->mask[i] = 0.0;
     444    }
    403445
    404446    return(newPoly);
    … …  
    408450{
    409451    int x = 0;
     452    int y = 0;
    410453    psDPolynomial2D *newPoly = NULL;
    411454
    … …  
    422465        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    423466    }
     467    for (x=0;x<nX;x++) {
     468        for (y=0;y<nY;y++) {
     469            newPoly->coeff[x][y] = 0.0;
     470            newPoly->coeffErr[x][y] = 0.0;
     471            newPoly->mask[x][y] = 0;
     472        }
     473    }
    424474
    425475    return(newPoly);
    … …  
    430480    int x = 0;
    431481    int y = 0;
     482    int z = 0;
    432483    psDPolynomial3D *newPoly = NULL;
    433484
    … …  
    450501        }
    451502    }
     503    for (x=0;x<nX;x++) {
     504        for (y=0;y<nY;y++) {
     505            for (z=0;z<nZ;z++) {
     506                newPoly->coeff[x][y][z] = 0.0;
     507                newPoly->coeffErr[x][y][z] = 0.0;
     508                newPoly->mask[x][y][z] = 0;
     509            }
     510        }
     511    }
    452512
    453513    return(newPoly);
    … …  
    459519    int x = 0;
    460520    int y = 0;
     521    int z = 0;
    461522    psDPolynomial4D *newPoly = NULL;
    462523
    … …  
    485546        }
    486547    }
     548    for (w=0;w<nW;w++) {
     549        for (x=0;x<nX;x++) {
     550            for (y=0;y<nY;y++) {
     551                for (z=0;z<nZ;z++) {
     552                    newPoly->coeff[w][x][y][z] = 0.0;
     553                    newPoly->coeffErr[w][x][y][z] = 0.0;
     554                    newPoly->mask[w][x][y][z] = 0;
     555                }
     556            }
     557        }
     558    }
    487559
    488560    return(newPoly);

  • trunk/psLib/src/math/psPolynomial.h

    r974 r1020  
     1/** @file psFunctions.h
     2 *  \brief Standard Mathematical Functions.
     3 *  \ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author George Gusciora, MHPCC
     13 *
     14 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2004-06-14 19:32:42 $
     16 *
     17 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     18 */
     19
    120#if !defined(PS_FUNCTIONS_H)
    221#define PS_FUNCTIONS_H
    322
    4 /** \file psFunctions.h
    5  *  \brief Standard Mathematical Functions.
    6  *  \ingroup Stats
    7  */
    823
    924/** \addtogroup Stats
    … …  
    2338
    2439
    25 psVector *psGaussianDev(float mean,
    26                         float sigma,
    27                         int Npts);
    28 
    29 /*****************************************************************************/
     40/** Produce a vector of random numbers from a Gaussian distribution with
     41    the specified mean and sigma */
     42psVector *psGaussianDev(float mean,    ///< The mean of the Gaussian
     43                        float sigma,   ///< The sigma of the Gaussian
     44                        int Npts);     ///< The size of the vector
     45
     46
     47
     48
    3049
    3150/** One-dimensional polynomial */
    3251typedef struct
    3352{
    34     int n;    ///< Number of terms
    35     float *coeff;   ///< Coefficients
    36     float *coeffErr;   ///< Error in coefficients
    37     char *mask;    ///< Coefficient mask
     53    int n;           ///< Number of terms
     54    float *coeff;    ///< Coefficients
     55    float *coeffErr; ///< Error in coefficients
     56    char *mask;      ///< Coefficient mask
    3857}
    3958psPolynomial1D;

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

    r889 r1020  
    44 *         routines.
    55 *
    6  *  @author George Gusciora, MHPCC
     6 *  This file will hold the functions for allocated, freeing, and evaluating
     7 *  polynomials.  It also contains a Gaussian functions.
    78 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-06-07 00:32:53 $
     9 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-06-14 19:32:42 $
    1011 *
    1112 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    103104psPolynomial1D *psPolynomial1DAlloc(int n)
    104105{
     106    int i = 0;
    105107    psPolynomial1D *newPoly = NULL;
    106108
    … …  
    110112    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
    111113    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     114    for (i=0;i<n;i++) {
     115        newPoly->coeff[i] = 0.0;
     116        newPoly->coeffErr[i] = 0.0;
     117        newPoly->mask[i] = 0.0;
     118    }
    112119
    113120    return(newPoly);
    … …  
    117124{
    118125    int x = 0;
     126    int y = 0;
    119127    psPolynomial2D *newPoly = NULL;
    120128
    … …  
    131139        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    132140    }
    133 
     141    for (x=0;x<nX;x++) {
     142        for (y=0;y<nY;y++) {
     143            newPoly->coeff[x][y] = 0.0;
     144            newPoly->coeffErr[x][y] = 0.0;
     145            newPoly->mask[x][y] = 0;
     146        }
     147    }
    134148
    135149    return(newPoly);
    … …  
    140154    int x = 0;
    141155    int y = 0;
     156    int z = 0;
    142157    psPolynomial3D *newPoly = NULL;
    143158
    … …  
    160175        }
    161176    }
     177    for (x=0;x<nX;x++) {
     178        for (y=0;y<nY;y++) {
     179            for (z=0;z<nZ;z++) {
     180                newPoly->coeff[x][y][z] = 0.0;
     181                newPoly->coeffErr[x][y][z] = 0.0;
     182                newPoly->mask[x][y][z] = 0;
     183            }
     184        }
     185    }
    162186
    163187    return(newPoly);
    … …  
    169193    int x = 0;
    170194    int y = 0;
     195    int z = 0;
    171196    psPolynomial4D *newPoly = NULL;
    172197
    … …  
    195220        }
    196221    }
     222    for (w=0;w<nW;w++) {
     223        for (x=0;x<nX;x++) {
     224            for (y=0;y<nY;y++) {
     225                for (z=0;z<nZ;z++) {
     226                    newPoly->coeff[w][x][y][z] = 0.0;
     227                    newPoly->coeffErr[w][x][y][z] = 0.0;
     228                    newPoly->mask[w][x][y][z] = 0;
     229                }
     230            }
     231        }
     232    }
    197233
    198234    return(newPoly);
    … …  
    394430psDPolynomial1D *psDPolynomial1DAlloc(int n)
    395431{
     432    int i = 0;
    396433    psDPolynomial1D *newPoly = NULL;
    397434
    … …  
    401438    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
    402439    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
     440    for (i=0;i<n;i++) {
     441        newPoly->coeff[i] = 0.0;
     442        newPoly->coeffErr[i] = 0.0;
     443        newPoly->mask[i] = 0.0;
     444    }
    403445
    404446    return(newPoly);
    … …  
    408450{
    409451    int x = 0;
     452    int y = 0;
    410453    psDPolynomial2D *newPoly = NULL;
    411454
    … …  
    422465        newPoly->mask[x]     = (char *)  psAlloc(nY * sizeof(char));
    423466    }
     467    for (x=0;x<nX;x++) {
     468        for (y=0;y<nY;y++) {
     469            newPoly->coeff[x][y] = 0.0;
     470            newPoly->coeffErr[x][y] = 0.0;
     471            newPoly->mask[x][y] = 0;
     472        }
     473    }
    424474
    425475    return(newPoly);
    … …  
    430480    int x = 0;
    431481    int y = 0;
     482    int z = 0;
    432483    psDPolynomial3D *newPoly = NULL;
    433484
    … …  
    450501        }
    451502    }
     503    for (x=0;x<nX;x++) {
     504        for (y=0;y<nY;y++) {
     505            for (z=0;z<nZ;z++) {
     506                newPoly->coeff[x][y][z] = 0.0;
     507                newPoly->coeffErr[x][y][z] = 0.0;
     508                newPoly->mask[x][y][z] = 0;
     509            }
     510        }
     511    }
    452512
    453513    return(newPoly);
    … …  
    459519    int x = 0;
    460520    int y = 0;
     521    int z = 0;
    461522    psDPolynomial4D *newPoly = NULL;
    462523
    … …  
    485546        }
    486547    }
     548    for (w=0;w<nW;w++) {
     549        for (x=0;x<nX;x++) {
     550            for (y=0;y<nY;y++) {
     551                for (z=0;z<nZ;z++) {
     552                    newPoly->coeff[w][x][y][z] = 0.0;
     553                    newPoly->coeffErr[w][x][y][z] = 0.0;
     554                    newPoly->mask[w][x][y][z] = 0;
     555                }
     556            }
     557        }
     558    }
    487559
    488560    return(newPoly);

  • trunk/psLib/src/math/psSpline.h

    r974 r1020  
     1/** @file psFunctions.h
     2 *  \brief Standard Mathematical Functions.
     3 *  \ingroup Stats
     4 *
     5 *  This file will hold the prototypes for procedures which allocate, free,
     6 *  and evaluate various polynomials.  Those polynomial structures are also
     7 *  defined here.
     8 *
     9 *  @ingroup Stats
     10 *
     11 *  @author Someone at IfA
     12 *  @author George Gusciora, MHPCC
     13 *
     14 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2004-06-14 19:32:42 $
     16 *
     17 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     18 */
     19
    120#if !defined(PS_FUNCTIONS_H)
    221#define PS_FUNCTIONS_H
    322
    4 /** \file psFunctions.h
    5  *  \brief Standard Mathematical Functions.
    6  *  \ingroup Stats
    7  */
    823
    924/** \addtogroup Stats
    … …  
    2338
    2439
    25 psVector *psGaussianDev(float mean,
    26                         float sigma,
    27                         int Npts);
    28 
    29 /*****************************************************************************/
     40/** Produce a vector of random numbers from a Gaussian distribution with
     41    the specified mean and sigma */
     42psVector *psGaussianDev(float mean,    ///< The mean of the Gaussian
     43                        float sigma,   ///< The sigma of the Gaussian
     44                        int Npts);     ///< The size of the vector
     45
     46
     47
     48
    3049
    3150/** One-dimensional polynomial */
    3251typedef struct
    3352{
    34     int n;    ///< Number of terms
    35     float *coeff;   ///< Coefficients
    36     float *coeffErr;   ///< Error in coefficients
    37     char *mask;    ///< Coefficient mask
     53    int n;           ///< Number of terms
     54    float *coeff;    ///< Coefficients
     55    float *coeffErr; ///< Error in coefficients
     56    char *mask;      ///< Coefficient mask
    3857}
    3958psPolynomial1D;

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

    r978 r1020  
    1515
    1616#include "float.h"
    17 #define ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
     17#define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
    1818// will use robust statistical methods.
    19 #define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar
    20 // smoothing.
    21 #define PI 3.141592653
     19#define GAUSS_WIDTH 20                  // The width of the Gaussian or boxcar smoothing.
    2220#define CLIPPED_NUM_ITER_LB 1
    2321#define CLIPPED_NUM_ITER_UB 10
    … …  
    2624#define true 1
    2725#define false 0
    28 #define MYMAXFLOAT 1e99
     26#define MYMAXFLOAT HUGE
    2927
    3028#ifndef DOXYGEN
    … …  
    4846    newStruct->sampleUQ = NAN;
    4947    newStruct->sampleLQ = NAN;
    50     newStruct->sampleLimit = ROBUST_SIZE_THRESHOLD;
     48    newStruct->sampleLimit = DEFAULT_ROBUST_SIZE_THRESHOLD;
    5149    newStruct->robustMean = NAN;
    5250    newStruct->robustMedian = NAN;
    … …  
    7876
    7977/******************************************************************************
    80     HISTOGRAM FUNCTIONS
     78psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure
     79with the specifed upper and lower limits, and the specifed number of bins.
     80This routine will also set the bounds for each of the bins.
     81 
     82Input:
     83    lower
     84    upper
     85    n
     86Returns:
     87    The histogram structure
    8188 *****************************************************************************/
    8289psHistogram *psHistogramAlloc(float lower,
    … …  
    8491                              int n)
    8592{
    86     int i = 0;
    87     psHistogram *newHist = NULL;
    88     float binSize = 0.0;
    89 
     93    int i = 0;                          // Loop index variable
     94    psHistogram *newHist = NULL;        // The new histogram structure
     95    float binSize = 0.0;                // The histogram bin size
     96
     97    // GUS: Verify that this is the correct action.
    9098    if (n == 0) {
    91         psAbort(__func__, "psHistogram requested with 0 bins");
     99        return(NULL);
    92100    }
    93101
    94102    // GUS: Verify that this is the correct action.
    95103    if (lower > upper) {
    96         //        psAbort(__func__, "psHistogram requested with lower>upper");
    97104        return(NULL);
    98105    }
    99106
     107    // Allocate memory for the new histogram structure.  If there are N
     108    // bins, then there are N+1 bounds to those bins.
    100109    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    101110    newHist->bounds = psVectorAlloc(n+1, PS_TYPE_F32);
    102111    newHist->bounds->n = newHist->bounds->nalloc;
    103112
     113    // Calculate the bounds for each bin.
    104114    binSize = (upper - lower) / (float) n;
    105115    for (i=0;i<n+1;i++) {
    106116        newHist->bounds->data.F32[i] = lower + (binSize * (float) i);
    107117    }
     118
     119    // Allocate the bins, and initialize them to zero.
    108120    newHist->nums =  psVectorAlloc(n, PS_TYPE_S32);
    109121    newHist->nums->n = newHist->nums->nalloc;
    … …  
    111123        newHist->nums->data.S32[i] = 0;
    112124    }
     125
     126    // Initialize the other members.
    113127    newHist->minNum = 0;
    114128    newHist->maxNum = 0;
    … …  
    118132}
    119133
    120 // When this is called, the number of data elements in bounds
    121 // is n.  Therefore, the number of bins is n-1.
     134/******************************************************************************
     135psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure
     136with the specifed bounds.  The number of elements in the bounds vector is n.
     137Therefore, the number of bins is n-1.
     138Input:
     139    bounds
     140Returns:
     141    The histogram structure
     142 *****************************************************************************/
    122143psHistogram *psHistogramAllocGeneric(const psVector *restrict bounds)
    123144{
    124     psHistogram *newHist = NULL;
    125     int i;
     145    psHistogram *newHist = NULL;        // The new histogram structure
     146    int i;                              // Loop index variable
    126147
    127148    // GUS: Verify that this is the correct action.
    … …  
    131152    }
    132153
     154    // GUS: Verify that this is the correct action.
     155    if (bounds->n <= 1) {
     156        //        psAbort(__func__, "psHistogram requested with NULL bounds");
     157        return(NULL);
     158    }
     159
     160    // Allocate memory for the new histogram structure.
    133161    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    134162    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
    … …  
    137165        newHist->bounds->data.F32[i] = bounds->data.F32[i];
    138166    }
     167
     168    // Allocate the bins, and initialize them to zero.  If there are N bounds,
     169    // then there are N-1 bins.
    139170    newHist->nums = psVectorAlloc((bounds->n)-1, PS_TYPE_S32);
    140171    newHist->nums->n = newHist->nums->nalloc;
    … …  
    143174    }
    144175
     176    // Initialize the other members.
    145177    newHist->minNum = 0;
    146178    newHist->maxNum = 0;
    … …  
    159191
    160192/*****************************************************************************
    161     NOTE: Can we assume fixed size bins in the myHist data structure?
     193psHistogramVector(out, in, mask, maskVal): this procedure takes as input a
     194preallocated and initialized histogram structure.  It fills the bins in that
     195histogram structure in accordance with the input data "in" and the, possibly
     196NULL, mask vector.
    162197 
    163     NOTE: This procedure assumes that the following has been already set
    164         myHist->lower
    165         myHist->upper
    166  myHist->minVal
    167  myHist->maxVal
    168     This procedure sets
    169         myHist->nums
    170         myHist->minNum
    171         myHist->maxNum
     198Inputs:
     199    out
     200    in
     201    mask
     202    maskVal
     203Returns:
     204    The histogram structure "out".
    172205 *****************************************************************************/
    173206psHistogram *psHistogramVector(psHistogram *out,
    174                                psVector *in,
    175                                psVector *mask,
     207                               const psVector *restrict in,
     208                               const psVector *restrict mask,
    176209                               int maskVal)
    177210{
    178     int i = 0;
    179     int j = 0;
    180     float binSize = 0.0;
    181     int binNum = 0;
    182     int numBins = 0;
     211    int i = 0;                                  // Loop index variable
     212    int j = 0;                                  // Loop index variable
     213    float binSize = 0.0;                        // Histogram bin size
     214    int binNum = 0;                             // A temporary bin number
     215    int numBins = 0;                            // The total number of bins
    183216
    184217    // GUS: Verify that this is the correct action.
    185218    if (out == NULL) {
    186         //        psAbort(__func__, "psHistogram requested with NULL struct");
    187219        return(NULL);
    188220    }
    … …  
    190222    // GUS: Verify that this is the correct action.
    191223    if (in == NULL) {
    192         //        psAbort(__func__, "psHistogram requested with NULL struct");
    193224        return(out);
    194225    }
    195226
     227    if (in->type.type != PS_TYPE_F32) {
     228        psAbort(__func__,
     229                "Only data type PS_TYPE_F32 is currently supported (0x%x).",
     230                in->type.type);
     231    }
     232
     233    if (mask != NULL) {
     234        if (in->n != mask->n) {
     235            psAbort(__func__,
     236                    "Vector data and vector mask are of different sizes.");
     237        }
     238        if (mask->type.type != PS_TYPE_U8) {
     239            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
     240        }
     241    }
     242    // GUS: determine the correct action for a variety of other cases:
     243    // in vector has 0 elements, and histogram structure has zero bins.
     244
    196245    numBins = out->nums->n;
    197246
    198     if (mask != NULL) {
    199         for (i=0;i<in->n;i++) {
    200             // Check if this pixel is masked, and if so, skip it.
    201             if (!(mask->data.U8[i] & maskVal)) {
    202                 // Check if this pixel is below the minimum value, and if so
    203                 // count it, then skip it.
    204                 if (in->data.F32[i] < out->bounds->data.F32[0]) {
    205                     out->minNum++;
    206 
    207                     // Check if this pixel is above the maximum value, and if so
    208                     // count it, then skip it.
    209                 } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    210                     out->maxNum++;
    211                 } else {
    212                     // If this is a uniform histogram, determining the correct
    213                     // number is trivial.
    214                     if (out->uniform == true) {
    215                         binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    216 
    217                         binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    218                                         binSize);
    219                         (out->nums->data.S32[binNum])++;
    220                         // If this is a non-uniform histogram, determining the correct
    221                         // bin number requires a bit more work.
    222                     } else {
    223                         // GUS: This is slow.  Put a smarter algorithm here to
    224                         // find the correct bin number (bin search, probably)
    225                         for (j=0;j<(out->bounds->n)-1;j++) {
    226                             if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
    227                                     (in->data.F32[i] <= out->bounds->data.S32[j+1])) {
    228                                 (out->nums->data.S32[j])++;
    229                             }
    230                         }
    231                     }
    232                 }
    233             }
    234         }
    235     } else {
    236         for (i=0;i<in->n;i++) {
     247    for (i=0;i<in->n;i++) {
     248        // Check if this pixel is masked, and if so, skip it.
     249        if ((mask == NULL) ||
     250                ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
     251            // Check if this pixel is below the minimum value, and if so
     252            // count it, then skip it.
    237253            if (in->data.F32[i] < out->bounds->data.F32[0]) {
    238                 // Check if this pixel is below the minimum value, and if so
    239                 // count it, then skip it.
    240254                out->minNum++;
    241             } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    242255                // Check if this pixel is above the maximum value, and if so
    243256                // count it, then skip it.
     257            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
    244258                out->maxNum++;
    245259            } else {
    … …  
    251265                    binNum = (int) ((in->data.F32[i] - out->bounds->data.F32[0]) /
    252266                                    binSize);
    253                     (out->nums->data.S32[binNum])++;
    254 
    255267                    // If this is a non-uniform histogram, determining the correct
    256268                    // bin number requires a bit more work.
    … …  
    267279            }
    268280        }
    269 
    270281    }
    271282    return(out);
    272283}
    273284
    274 void p_printVector(psVector *myVector,
    275                    psVector *maskVector,
    276                    unsigned int maskVal,
    277                    psStats *stats)
    278 {
    279     int i = 0;
     285/*****************************************************************************
     286p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
     287function that simply prints a vector to STDOUT.  Used primarily for
     288debugging.
     289 *****************************************************************************/
     290
     291void p_psVectorPrint(psVector *myVector,
     292                     psVector *maskVector,
     293                     unsigned int maskVal,
     294                     psStats *stats)
     295{
     296    int i = 0;                                  // Loop index variable.
    280297
    281298    for (i=0;i<myVector->n;i++) {
    … …  
    291308 ******************************************************************************
    292309 ******************************************************************************
    293     MISC STATISTICAL FUNCTIONS
     310    MISC PRIVATE STATISTICAL FUNCTIONS
     311 
     312    NOTE: it is assumed that any call to these statistical functions will
     313    have been preceded by a call to the psVectorStats() function.  Various
     314    sanity tests will only be performed in psVectorStats().
     315        Is the mask vector the same length as the data vector?
     316        Is the mask vector of type PS_TYPE_U8?
     317        Is the stats data structure NULL?
     318        Is the in data structure NULL?
     319        Is the in data structure of type PS_TYPE_F32?
    294320 ******************************************************************************
    295321 ******************************************************************************
    … …  
    298324
    299325/******************************************************************************
    300     ASSUMPTION: the mean is always calculated exactly.  Robust means are never
    301     calculated in this routine.
     326p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
     327mean of the input vector.
     328Inputs
     329Returns
     330    NULL
     331ASSUMPTION: the mean is always calculated exactly.  Robust means are never
     332calculated in this routine.
    302333 *****************************************************************************/
    303334void p_psVectorSampleMean(const psVector *restrict myVector,
    … …  
    306337                          psStats *stats)
    307338{
    308     int i = 0;
    309     float mean = 0.0;
    310     int count = 0;
    311     float rangeMin = 0.0;
    312     float rangeMax = 0.0;
    313 
     339    int i = 0;                                  // Loop index variable
     340    float mean = 0.0;                           // The mean
     341    int count = 0;                              // # of points in this mean?
     342    float rangeMin = 0.0;                       // Exclude data below this
     343    float rangeMax = 0.0;                       // Exclude date above this
     344
     345    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
     346    // loop.
    314347    if (stats->options & PS_STAT_USE_RANGE) {
    315348        rangeMin = stats->min;
    … …  
    317350        if (maskVector != NULL) {
    318351            for (i=0;i<myVector->n;i++) {
     352                // Check if the data is with the specified range
    319353                if (!(maskVal & maskVector->data.U8[i]) &&
    320354                        (rangeMin <= myVector->data.F32[i]) &&
    … …  
    355389}
    356390
     391/******************************************************************************
     392p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     393max of the input vector.
     394Inputs
     395    myVector
     396    maskVector
     397    maskVal
     398    stats
     399Returns
     400    NULL
     401 *****************************************************************************/
    357402void p_psVectorMax(const psVector *restrict myVector,
    358403                   const psVector *restrict maskVector,
    … …  
    360405                   psStats *stats)
    361406{
    362     int i = 0;
    363     float max = -MYMAXFLOAT;
    364     float rangeMin = 0.0;
    365     float rangeMax = 0.0;
    366 
     407    int i = 0;                                  // Loop index variable
     408    float max = -MYMAXFLOAT;                    // The calculated maximum
     409    float rangeMin = 0.0;                       // Exclude data below this
     410    float rangeMax = 0.0;                       // Exclude date above this
     411
     412    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    367413    if (stats->options & PS_STAT_USE_RANGE) {
    368414        rangeMin = stats->min;
    … …  
    408454}
    409455
     456/******************************************************************************
     457p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
     458minimum of the input vector.
     459Inputs
     460    myVector
     461    maskVector
     462    maskVal
     463    stats
     464Returns
     465    NULL
     466 *****************************************************************************/
    410467void p_psVectorMin(const psVector *restrict myVector,
    411468                   const psVector *restrict maskVector,
    … …  
    413470                   psStats *stats)
    414471{
    415     int i = 0;
    416     float min = MYMAXFLOAT;
    417     float rangeMin = 0.0;
    418     float rangeMax = 0.0;
     472    int i = 0;                                  // Loop index variable
     473    float min = MYMAXFLOAT;                     // The calculated maximum
     474    float rangeMin = 0.0;                       // Exclude data below this
     475    float rangeMax = 0.0;                       // Exclude date above this
    419476
    420477    if (stats->options & PS_STAT_USE_RANGE) {
    … …  
    462519
    463520/******************************************************************************
    464     This routine calculates the number of non-masked pixels in the vector.
    465     that fall within the min/max range, if given.
     521p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
     522number of non-masked pixels in the vector.  that fall within the min/max
     523range, if given.
     524Inputs
     525    myVector
     526    maskVector
     527    maskVal
     528    stats
     529Returns
     530    NULL
    466531 *****************************************************************************/
    467532int p_psVectorNValues(const psVector *restrict myVector,
    … …  
    470535                      psStats *stats)
    471536{
    472     int i = 0;
    473     int numData = 0;
    474     float rangeMin = 0.0;
    475     float rangeMax = 0.0;
     537    int i = 0;                                  // Loop index variable
     538    int numData = 0;                            // The number of data points
     539    float rangeMin = 0.0;                       // Exclude data below this
     540    float rangeMax = 0.0;                       // Exclude date above this
    476541
    477542    if (stats->options & PS_STAT_USE_RANGE) {
    … …  
    512577
    513578
     579/******************************************************************************
     580p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
     581median of the input vector.
     582Inputs
     583    myVector
     584    maskVector
     585    maskVal
     586    stats
     587Returns
     588    NULL
     589 *****************************************************************************/
    514590void p_psVectorSampleMedian(const psVector *restrict myVector,
    515591                            const psVector *restrict maskVector,
    … …  
    517593                            psStats *stats)
    518594{
    519     psVector *unsortedVector = NULL;
    520     psVector *sortedVector = NULL;
    521     int count = 0;
    522     int i = 0;
    523     int nValues = 0;
    524     float rangeMin = 0.0;
    525     float rangeMax = 0.0;
    526     psStats *stats2 = NULL;
     595    psVector *unsortedVector = NULL;            // Temporary vector
     596    psVector *sortedVector = NULL;              // Temporary vector
     597    int i = 0;                                  // Loop index variable
     598    int count = 0;                              // # of points in this mean?
     599    int nValues = 0;                            // # of points in vector
     600    float rangeMin = 0.0;                       // Exclude data below this
     601    float rangeMax = 0.0;                       // Exclude date above this
     602    psStats *stats2 = NULL;                     // Temporary stats structure
    527603
    528604
    … …  
    531607
    532608    if (myVector->n > stats->sampleLimit) {
     609        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     610
    533611        // Calculate the robust quartiles.
    534612        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    … …  
    548626
    549627    // Determine how many data points fit inside this min/max range
    550     // and are not maxed, IF the maskVector is not NULL>
     628    // and are not masked, IF the maskVector is not NULL>
    551629    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    552630
    … …  
    620698                                  float sigma)
    621699{
    622     int i = 0;
    623     int j = 0;
    624     float tmpf = 0.0;
    625     float gaussianCoefs[1 + (2 * GAUSS_WIDTH)];
    626     // The coefficients used in the histogram
    627     // smoothing calculation.
     700    int i = 0;                                  // Loop index variable
     701    int j = 0;                                  // Loop index variable
     702    float tmpf = 0.0;                           // Temporary variable
     703    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
    628704
    629705    for(i=0;i<(1 + (2 * GAUSS_WIDTH));i++) {
    … …  
    631707            // If sigma does not equal zero, then we use Gaussian smoothing.
    632708            #ifdef  DARWIN
    633             tmpf = (float) sqrt(2.0f * PI * sigma * sigma);
     709            tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma);
    634710            #else
    635711
    636             tmpf = sqrtf(2.0f * PI * sigma * sigma);
     712            tmpf = sqrtf(2.0f * M_PI * sigma * sigma);
    637713            #endif
    638714
    639715            gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) *
    640716                                             ((float) (i-GAUSS_WIDTH))) /
    641                                             (2.0f * sigma * sigma)) /
    642                                tmpf;
     717                                            (2.0f * sigma * sigma)) / tmpf;
    643718        } else {
    644719            /* If sigma equals zero (all pixels have the same value)
    645              * the above code will divide by zero.  Therefore, we instead
    646              * use boxcar smoothing.
     720             * the above code will divide by zero.  Therefore, we don't need
     721             * to smooth the data.
    647722             */
    648             gaussianCoefs[i] = 1.0f / (1.0f + (2.0f * (float) GAUSS_WIDTH));
     723            return;
    649724        }
    650725    }
    … …  
    662737
    663738/******************************************************************************
    664     p_psVectorSampleQuartiles()
    665  This procedure calculates the upper and/or lower quartiles of the
    666  data set.
    667  
     739p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
     740the upper and/or lower quartiles of the input vector.
     741Inputs
     742    myVector
     743    maskVector
     744    maskVal
     745    stats
     746Returns
     747    NULL
    668748 *****************************************************************************/
    669749void p_psVectorSampleQuartiles(const psVector *restrict myVector,
    … …  
    672752                               psStats *stats)
    673753{
    674     psVector *unsortedVector = NULL;
    675     psVector *sortedVector = NULL;
    676     int count = 0;
    677     int ind = 0;
    678     int i = 0;
    679     int nValues = 0;
    680     float rangeMin = 0.0;
    681     float rangeMax = 0.0;
     754    psVector *unsortedVector = NULL;            // Temporary vector
     755    psVector *sortedVector = NULL;              // Temporary vector
     756    int i = 0;                                  // Loop index variable
     757    int count = 0;                              // # of points in this mean?
     758    int nValues = 0;                            // # data points
     759    float rangeMin = 0.0;                       // Exclude data below this
     760    float rangeMax = 0.0;                       // Exclude date above this
    682761
    683762    // Determine if the number of data points exceed a threshold which will
    684763    // cause to generate robust stats, as opposed to exact stats.
    685     /* GUS: When IfA provides the algorithm, insert is here.
    686         if (myVector->n > stats->sampleLimit) {
    687             psStats *stats2 = NULL;
    688             // Calculate the robust quartiles.
    689             stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
    690             p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    691      
    692             // Store the robust quartiles into the sample quartile members.
    693             stats->sampleUQ = stats2->robustUQ;
    694             stats->sampleLQ = stats2->robustLQ;
    695      
    696             // Free temporary data buffers.
    697             psStatsFree(stats2);
    698      
    699             // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    700             stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    701      
    702             return;
    703         }
    704     */
     764    if (myVector->n > stats->sampleLimit) {
     765        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     766        psStats *stats2 = NULL;
     767        // Calculate the robust quartiles.
     768        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
     769        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     770
     771        // Store the robust quartiles into the sample quartile members.
     772        stats->sampleUQ = stats2->robustUQ;
     773        stats->sampleLQ = stats2->robustLQ;
     774
     775        // Free temporary data buffers.
     776        psStatsFree(stats2);
     777
     778        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     779        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     780
     781        return;
     782    }
    705783
    706784    // Determine how many data points fit inside this min/max range
    … …  
    757835
    758836    // Calculate the quartile points exactly.
    759     ind = 3 * (nValues / 4);
    760     stats->sampleUQ = sortedVector->data.F32[ind];
    761     ind = (nValues / 4);
    762     stats->sampleLQ = sortedVector->data.F32[ind];
     837    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     838    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    763839
    764840    // Free the temporary data structures.
    … …  
    770846
    771847/******************************************************************************
    772     p_psVectorRobustStats(): this procedure calculates a variety of robust
    773     stat measures:
     848p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
     849calculates a variety of robust stat measures:
    774850        PS_STAT_ROBUST_MEAN
    775851        PS_STAT_ROBUST_MEDIAN
    … …  
    784860    measures were calcualted in separate functiosn, then much of the initial
    785861    processing would be duplicated.
     862Inputs
     863    myVector
     864    maskVector
     865    maskVal
     866    stats
     867Returns
     868    NULL
    786869 *****************************************************************************/
    787870void p_psVectorRobustStats(const psVector *restrict myVector,
    … …  
    791874{
    792875    psHistogram *robustHistogram = NULL;
    793     float binSize = 0.0;
     876    float binSize = 0.0;                        // Size of the histogram bins
    794877    float sigmaE = 0.0;
    795     int LQBinNum = -1;
    796     int UQBinNum = -1;
    797     int i = 0;    // Loop index variable.
     878    int LQBinNum = -1;                          // Bin num for lower quartile
     879    int UQBinNum = -1;                          // Bin num for upper quartile
     880    int i = 0;                                  // Loop index variable
    798881    int maxBinNum = 0;
    799882    int maxBinCount = 0;
    800883    float dL = 0.0;
     884    int numBins = 0;
    801885
    802886    // NOTE: The SDRS states that the sample quartiles must be used to
    … …  
    827911    }
    828912
    829     // Create the histogram structure.
     913    // Create the histogram structure (yes, 2 is necessary, not 1).
     914    numBins = 2 + (int) ((stats->max - stats->min) / binSize);
     915
    830916    robustHistogram = psHistogramAlloc(stats->min,
    831917                                       stats->max,
    832                                        1 + (int) ((stats->max - stats->min) / binSize));
    833 
    834 
     918                                       numBins);
    835919    // Populate the histogram array.
    836     // GUS: fix this
    837     //    robustHistogram = psHistogramVector(robustHistogram, myVector);
    838     //
     920    psHistogramVector(robustHistogram, myVector, maskVector, maskVal);
    839921
    840922    // Smooth the histogram.
    … …  
    868950
    869951    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    870     // What algorithm should I use to do this?
    871 
    872     // The resulting fit parameters are the robust mean, mean_r, and sigma
    873     // What is the mean_r?
     952    // GUS: This step is dependent on the functions in psMinimize.c being
     953    // implemented.  Currently, they are not.
    874954
    875955    if (stats->options & PS_STAT_ROBUST_MEAN) {
    … …  
    895975}
    896976
    897 /*****************************************************************************
    898     NOTE: This function assumes that p_psVectorMean() has already been called
    899     and the correct value is stored in stats->sampleMean.
     977/******************************************************************************
     978p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
     979stdev of the input vector.
     980Inputs
     981    myVector
     982    maskVector
     983    maskVal
     984    stats
     985Returns
     986    NULL
    900987 
    901     NOTE: the mean is always calculated exactly.  Robust means are never
    902     calculated in this routine.
     988NOTE: the mean is always calculated exactly.  Robust means are never
     989calculated in this routine.
    903990 *****************************************************************************/
    904991void p_psVectorSampleStdev(const psVector *restrict myVector,
    … …  
    907994                           psStats *stats)
    908995{
    909     int i            = 0;
    910     int countInt     = 0;
    911     float countFloat = 0.0;
    912     float       mean = 0.0;
    913     float       diff = 0.0;
    914     float sumSquares = 0.0;
    915     float sumDiffs = 0.0;
    916     float rangeMin = 0.0;
    917     float rangeMax = 0.0;
    918 
     996    int i = 0;                                  // Loop index variable
     997    int countInt = 0;                           // # of data points being used
     998    float countFloat = 0.0;                     // # of data points being used
     999    float mean = 0.0;                           // The mean
     1000    float diff = 0.0;                           // Used in calculating stdev
     1001    float sumSquares = 0.0;                     // temporary variable
     1002    float sumDiffs = 0.0;                       // temporary variable
     1003    float rangeMin = 0.0;                       // Exclude data below this
     1004    float rangeMax = 0.0;                       // Exclude date above this
     1005
     1006    // This procedure requires the mean.  If it has not been already
     1007    // calculated, then call p_psVectorSampleMean()
    9191008    if (0 != isnan(stats->sampleMean)) {
    9201009        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    … …  
    9791068}
    9801069
     1070/******************************************************************************
     1071p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
     1072clipped stats (mean or stdev) of the input vector.
     1073 
     1074Inputs
     1075    myVector
     1076    maskVector
     1077    maskVal
     1078    stats
     1079Returns
     1080    NULL
     1081 *****************************************************************************/
    9811082void p_psVectorClippedStats(const psVector *restrict myVector,
    9821083                            const psVector *restrict maskVector,
    … …  
    9841085                            psStats *stats)
    9851086{
    986     int i = 0;
    987     int j = 0;
    988     float clippedMean = 0.0;
    989     float clippedStdev = 0.0;
    990     psVector *tmpMask = NULL;
    991 
     1087    int i = 0;                                  // Loop index variable
     1088    int j = 0;                                  // Loop index variable
     1089    float clippedMean = 0.0;                    // self-explanatory
     1090    float clippedStdev = 0.0;                   // self-explanatory
     1091    float oldStanMean = 0.0;                    // Temporary variable
     1092    float oldStanStdev = 0.0;                   // Temporary variable
     1093    psVector *tmpMask = NULL;                   // Temporary vector
     1094
     1095    // Endure that stats->clipIter is within the proper range.
    9921096    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
    9931097            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
    … …  
    9961100    }
    9971101
     1102    // Endure that stats->clipSigma is within the proper range.
    9981103    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
    9991104            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
    … …  
    10021107    }
    10031108
    1004 
     1109    // We allocate a temporary mask vector since during the iterative
     1110    // steps that follow, we will be masking off additional data points.
     1111    // However, we do no want to modify the original mask vector.
    10051112    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    10061113    tmpMask->n = myVector->n;
    10071114
     1115    // If we were called with a mask vector, then initialize the temporary
     1116    // mask vector with those values.
    10081117    if (maskVector != NULL) {
    10091118        for (i=0;i<tmpMask->n;i++) {
    … …  
    10261135    clippedStdev = stats->sampleStdev;
    10271136
     1137    // Must save the old sampleMean and sampleStdev since the following code
     1138    // block overwrites them.
     1139    oldStanMean = stats->sampleMean;
     1140    oldStanStdev = stats->sampleStdev;
     1141
    10281142    // 5. Repeat N times:
    1029 
    10301143    for (i=0;i<stats->clipIter;i++) {
    1031         //        printf("p_psVectorClippedStats(): Iteration %d (%f %f)\n", i,
    1032         //                clippedMean, clippedStdev);
    10331144        for (j=0;j<myVector->n;j++) {
    10341145            // a) Exclude all values x_i for which |x_i - x| > K * stdev
    … …  
    10381149            }
    10391150            // b) compute new mean and stdev
    1040             // GUS: I should probably create a new struct here since the
    1041             // following calls will overwrite any old values in sampleMean.
    10421151            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
    10431152            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
    … …  
    10501159        }
    10511160    }
     1161    stats->sampleMean = oldStanMean;
     1162    stats->sampleStdev= oldStanStdev;
    10521163
    10531164    // 7. The last calcuated value of x is the cliped mean.
    … …  
    10661177
    10671178/******************************************************************************
     1179psVectorStats(myVector, maskVector, maskVal, stats): this is the public API
     1180function which calls the above private stats functions based on what bits
     1181were set in stats->options.
    10681182 
    1069     NOTE: The current strategy is to implement everything assuming that all
    1070     input data is of type PS_TYPE_F32.  Once the basic code is in place,
    1071     we will macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
     1183Inputs
     1184    myVector
     1185    maskVector
     1186    maskVal
     1187    stats
     1188Returns
     1189    The stats structure.
    10721190 
     1191NOTE: The current strategy is to implement everything assuming that all
     1192input data is of type PS_TYPE_F32.  Once the basic code is in place, we will
     1193macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
    10731194 *****************************************************************************/
    10741195psStats *psVectorStats(psStats *stats,
    … …  
    10771198                       unsigned int maskVal)
    10781199{
     1200    // GUS: Verify that this is the correct action.
    10791201    if (in == NULL) {
    10801202        return(stats);
    … …  
    10841206    }
    10851207
     1208    // Ensure that the data is of type PS_TYPE_F32.  Eventually, more data
     1209    // types will be implemented.
    10861210    if (in->type.type != PS_TYPE_F32) {
    10871211        psAbort(__func__,
    … …  
    10901214    }
    10911215
     1216    // Ensure that the mask vector is of the proper size and type.
    10921217    if (mask != NULL) {
    10931218        if (in->n != mask->n) {
    … …  
    11241249    }
    11251250
     1251    // Since the various robust stats quantities share much computation, they
     1252    // are grouped together in a single private function:
     1253    // p_psVectorRobustStats()
    11261254    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
    11271255            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
    … …  
    11471275    }
    11481276
    1149     //    OLD CODE: Should we check for an unknown option?
    1150     //    default:
    1151     //        psAbort(__func__, "Unknown options 0x%x.\n", stats->options);
    1152 
    11531277    return(stats);
    11541278}

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

    r974 r1020  
     1/** @file  psStats.h
     2 *  \brief basic statistical operations
     3 *  @ingroup Stats
     4 *
     5 *  This file will hold the definition of the histogram and stats data
     6 *  structures.  It also contains prototypes for procedures which operate
     7 *  on those data structures.
     8 *
     9 *  @author George Gusciora, MHPCC
     10 *
     11 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-06-14 19:33:09 $
     13 *
     14 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     15 */
    116#if !defined(PS_STATS_H)
    217#define PS_STATS_H
    3 /******************************************************************************
    4     This file will histogram/statistical functions and data structures.
    5  *****************************************************************************/
    6 /** \file psStats.h
    7  *  \brief basic statistical operations
    8  *  \ingroup Stats
    9  */
     18
    1019#include "psVector.h"
    1120
    … …  
    3645
    3746
    38 /** This is the generic statistics structure */
     47/** This is the generic statistics structure.  It contails the data members
     48    for the various statistic values.  It contains the options member to
     49    specifiy which statistics must be calculated. */
    3950typedef struct
    4051{
    … …  
    7485             );
    7586
    76 /** Constructor */
     87/** A constructor for the stats structure.*/
    7788psStats *psStatsAlloc(psStatsOptions options); ///< Statistics to measure
    7889
    79 /** Destructor */
     90/** A destructor for the stats structure.*/
    8091void psStatsFree(psStats *restrict stats); ///< Stats structure to destroy
    8192
    … …  
    8394    Histogram functions and data structures.
    8495 *****************************************************************************/
     96/** The basic histogram structure which contains bounds and bins. */
    8597typedef struct
    8698{
    … …  
    111123/** Calculate a histogram \ingroup MathGroup **/
    112124psHistogram *psHistogramVector (psHistogram *out,   ///< Histogram data
    113                                 psVector *in,       ///< Vector to analyse
    114                                 psVector *mask,     ///< Mask dat for input vector
     125                                const psVector *restrict in,       ///< Vector to analyse
     126                                const psVector *restrict mask,     ///< Mask dat for input vector
    115127                                int maskVal);       ///< Mask value
    116128
Note: See TracChangeset for help on using the changeset viewer.