IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 652


Ignore:
Timestamp:
May 12, 2004, 3:05:35 PM (22 years ago)
Author:
gusciora
Message:

Added the median calculation functions.

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r643 r652  
    33#include <string.h>
    44#include <stdarg.h>
    5 #include "pslib.h"
    6 /******************************************************************************
     5/******************************************************************************
     6//#include "pslib.h"
    77//#include "psMemory.h"
    88//#include "psTrace.h"
    … …  
    6464    psHistogram *newHist = NULL;
    6565    int numBins = 0;
     66    psType dummyFloatType;
     67    psType dummyIntType;
     68
     69    dummyFloatType.type = PS_TYPE_FLOAT;
     70    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     71    dummyIntType.type = PS_TYPE_FLOAT;
     72    dummyIntType.dimen = PS_DIMEN_VECTOR;
    6673
    6774    numBins = 1 + ((upper - lower) / size);
    6875    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    69     newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    70     newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    71     newHist->nums =  psVectorAlloc(numBins, PS_TYPE_INT);
     76    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
     77    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
     78    newHist->nums =  psVectorAlloc(dummyIntType, numBins);
    7279    newHist->minVal = lower;
    7380    newHist->maxVal = upper;
    … …  
    8693    int numBins = 0;
    8794    int i;
     95    psType dummyFloatType;
     96    psType dummyIntType;
     97
     98    dummyFloatType.type = PS_TYPE_FLOAT;
     99    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     100    dummyIntType.type = PS_TYPE_FLOAT;
     101    dummyIntType.dimen = PS_DIMEN_VECTOR;
    88102
    89103    if (lower->n != upper->n) {
    … …  
    92106    numBins = lower->n;
    93107    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    94     newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    95     newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
     108    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
     109    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
    96110    for (i=0;i<numBins;i++) {
    97         newHist->lower->vec.vecF[i] = lower->vec.vecF[i];
    98         newHist->upper->vec.vecF[i] = upper->vec.vecF[i];
    99     }
    100     newHist->nums = psVectorAlloc(numBins, PS_TYPE_INT);
     111        newHist->lower->vec.f[i] = lower->vec.f[i];
     112        newHist->upper->vec.f[i] = upper->vec.f[i];
     113    }
     114    newHist->nums = psVectorAlloc(dummyIntType, numBins);
    101115    newHist->minVal = minVal;
    102116    newHist->maxVal = maxVal;
    … …  
    135149    binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
    136150    for (i=0;i<myVector->n;i++) {
    137         if (myVector->vec.vecF[i] < myHist->minVal) {
     151        if (myVector->vec.f[i] < myHist->minVal) {
    138152            myHist->minNum++;
    139         } else if (myVector->vec.vecF[i] > myHist->maxVal) {
     153        } else if (myVector->vec.f[i] > myHist->maxVal) {
    140154            myHist->maxNum++;
    141155        } else {
    142             binNum = (int) ((myVector->vec.vecF[i] - myHist->minVal) / binSize);
    143             if ((myVector->vec.vecF[i] >= myHist->lower->vec.vecF[i]) &&
    144                     (myVector->vec.vecF[i] <= myHist->upper->vec.vecF[i])) {
    145                 myHist->nums->vec.vecI[i]++;
     156            binNum = (int) ((myVector->vec.f[i] - myHist->minVal) / binSize);
     157            if ((myVector->vec.f[i] >= myHist->lower->vec.f[i]) &&
     158                    (myVector->vec.f[i] <= myHist->upper->vec.f[i])) {
     159                myHist->nums->vec.i[i]++;
    146160            } else {
    147161                psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
    … …  
    164178    if (maskVector != NULL) {
    165179        for (i=0;i<myVector->n;i++) {
    166             if (!(maskVal & maskVector->vec.vecI[i])) {
    167                 mean+= maskVector->vec.vecI[i];
     180            if (!(maskVal & maskVector->vec.i[i])) {
     181                mean+= maskVector->vec.i[i];
    168182            }
    169183        }
    170184    } else {
    171185        for (i=0;i<myVector->n;i++) {
    172             mean+= maskVector->vec.vecI[i];
     186            mean+= maskVector->vec.i[i];
    173187        }
    174188    }
    … …  
    185199    if (maskVector != NULL) {
    186200        for (i=0;i<myVector->n;i++) {
    187             if (!(maskVal & maskVector->vec.vecI[i])) {
    188                 if (myVector->vec.vecF[i] > max) {
    189                     max = maskVector->vec.vecI[i];
     201            if (!(maskVal & maskVector->vec.i[i])) {
     202                if (myVector->vec.f[i] > max) {
     203                    max = maskVector->vec.i[i];
    190204                }
    191205            }
    … …  
    193207    } else {
    194208        for (i=0;i<myVector->n;i++) {
    195             if (myVector->vec.vecF[i] > max) {
    196                 max = maskVector->vec.vecI[i];
     209            if (myVector->vec.f[i] > max) {
     210                max = maskVector->vec.i[i];
    197211            }
    198212        }
    … …  
    210224    if (maskVector != NULL) {
    211225        for (i=0;i<myVector->n;i++) {
    212             if (!(maskVal & maskVector->vec.vecI[i])) {
    213                 if (myVector->vec.vecF[i] < min) {
    214                     min = maskVector->vec.vecI[i];
     226            if (!(maskVal & maskVector->vec.i[i])) {
     227                if (myVector->vec.f[i] < min) {
     228                    min = maskVector->vec.i[i];
    215229                }
    216230            }
    … …  
    218232    } else {
    219233        for (i=0;i<myVector->n;i++) {
    220             if (myVector->vec.vecF[i] < min) {
    221                 min = maskVector->vec.vecI[i];
     234            if (myVector->vec.f[i] < min) {
     235                min = maskVector->vec.i[i];
    222236            }
    223237        }
    … …  
    226240}
    227241
     242/******************************************************************************
     243 
     244 *****************************************************************************/
    228245int p_psArrayNValues(const psVector *restrict myVector, // float
    229246                     const psVector *restrict maskVector, // ints
    … …  
    235252    if (maskVector != NULL) {
    236253        for (i=0;i<myVector->n;i++) {
    237             if (!(maskVal & maskVector->vec.vecI[i])) {
     254            if (!(maskVal & maskVector->vec.i[i])) {
    238255                numData++;
    239256            }
    … …  
    243260    }
    244261    return(numData);
     262}
     263
     264
     265
     266#define MEDIAN_SIZE_THRESHOLD 10000
     267float p_psArraySampleMedian(const psVector *restrict myVector,
     268                            const psVector *restrict maskVector,
     269                            unsigned int maskVal)
     270{
     271    psVector *unsortedVector = NULL;
     272    psVector *sortedVector = NULL;
     273    psType dummyFloatType;
     274    int dataSize = 0;
     275    int count = 0;
     276    int i = 0;
     277    float median = 0.0;
     278
     279    dummyFloatType.type = PS_TYPE_FLOAT;
     280    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     281    dataSize = p_psArrayNValues(myVector, maskVector, maskVal);
     282
     283    if (dataSize < MEDIAN_SIZE_THRESHOLD) {
     284        unsortedVector = psVectorAlloc(dummyFloatType, dataSize);
     285        sortedVector   = psVectorAlloc(dummyFloatType, dataSize);
     286
     287        count = 0;
     288        if (maskVector != NULL) {
     289            for (i=0;i<myVector->n;i++) {
     290                if (!(maskVal & maskVector->vec.i[i])) {
     291                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     292                }
     293            }
     294            psSort(sortedVector, unsortedVector);
     295        } else {
     296            psSort(sortedVector, myVector);
     297        }
     298
     299        if (0 == (dataSize % 2)) {
     300            median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
     301                            sortedVector->vec.f[dataSize/2]);
     302        } else {
     303            median = sortedVector->vec.f[dataSize/2];
     304        }
     305    } else {
     306        // BROAD: Calculate the Robust Median
     307        // Determine the LQ of the distribution.
     308        // Determine the UQ of the distribution.
     309        // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
     310        // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
     311        // Find the bin with the peak value between LQ and UQ (the MODE)
     312        // dL = (UQ - LQ) / 8
     313        // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     314        // The resulting fit parameters are the robust mean, mean_r, and sigma
     315        psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
     316    }
     317
     318    psVectorFree(unsortedVector);
     319    psVectorFree(sortedVector);
     320    return(median);
     321}
     322
     323
     324int p_psArrayXXX(const psVector *restrict myVector,
     325                 const psVector *restrict maskVector,
     326                 unsigned int maskVal)
     327{
     328    printf("ERROR: Don't call me: p_psArrayXXX()\n");
     329    exit(1);
    245330}
    246331
    … …  
    265350    switch (stats->options) {
    266351    case PS_STAT_SAMPLE_MEAN:
    267         newStruct->max = p_psArraySampleMean(myVector, maskVector, maskVal);
     352        newStruct->sampleMean = p_psArraySampleMean(myVector, maskVector, maskVal);
    268353        break;
    269354    case PS_STAT_SAMPLE_MEDIAN:
    270         printf("Code me!\n");
     355        newStruct->sampleMedian = p_psArrayXXX(myVector, maskVector, maskVal);
    271356        break;
    272357    case PS_STAT_SAMPLE_STDEV:
    273         printf("Code me!\n");
     358        newStruct->sampleStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    274359        break;
    275360    case PS_STAT_SAMPLE_UQ:
    276         printf("Code me!\n");
     361        newStruct->sampleUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    277362        break;
    278363    case PS_STAT_SAMPLE_LQ:
    279         printf("Code me!\n");
     364        newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    280365        break;
    281366    case PS_STAT_ROBUST_MEAN:
    282         printf("Code me!\n");
     367        newStruct->robustMean = p_psArrayXXX(myVector, maskVector, maskVal);
    283368        break;
    284369    case PS_STAT_ROBUST_MEAN_NVALUES:
    285         printf("Code me!\n");
     370        newStruct->robustMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    286371        break;
    287372    case PS_STAT_ROBUST_MEDIAN:
    288         printf("Code me!\n");
     373        newStruct->robustMedian = p_psArrayXXX(myVector, maskVector, maskVal);
    289374        break;
    290375    case PS_STAT_ROBUST_MEDIAN_NVALUES:
    291         printf("Code me!\n");
     376        newStruct->robustMedianNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    292377        break;
    293378    case PS_STAT_ROBUST_MODE:
    294         printf("Code me!\n");
     379        newStruct->robustMode = p_psArrayXXX(myVector, maskVector, maskVal);
    295380        break;
    296381    case PS_STAT_ROBUST_MODE_NVALUES:
    297         printf("Code me!\n");
     382        newStruct->robustModeNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    298383        break;
    299384    case PS_STAT_ROBUST_STDEV:
    300         printf("Code me!\n");
     385        newStruct->robustStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    301386        break;
    302387    case PS_STAT_ROBUST_UQ:
    303         printf("Code me!\n");
     388        newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    304389        break;
    305390    case PS_STAT_ROBUST_LQ:
    306         printf("Code me!\n");
     391        newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    307392        break;
    308393    case PS_STAT_CLIPPED_MEAN:
    309         printf("Code me!\n");
     394        newStruct->clippedMean = p_psArrayXXX(myVector, maskVector, maskVal);
    310395        break;
    311396    case PS_STAT_CLIPPED_MEAN_NVALUES:
    312         printf("Code me!\n");
     397        newStruct->clippedMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    313398        break;
    314399    case PS_STAT_CLIPPED_MEAN_NSIGMA:
    315         printf("Code me!\n");
     400        newStruct->clipSigma = p_psArrayXXX(myVector, maskVector, maskVal);
    316401        break;
    317402    case PS_STAT_CLIPPED_STDEV:
    318         printf("Code me!\n");
     403        newStruct->clippedStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    319404        break;
    320405    case PS_STAT_MAX:

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

    r644 r652  
    66 *  \ingroup MathGroup
    77 */
     8#include "psLib.h"
     9#include "psVector.h"
    810
    9 /******************************************************************************
    10     The following typedefs and functions belong in the psArray.h and psArray.c
    11     files.  However, those files are not available at this time, due to the
    12     changing definition of the psLib data/functions.  For temporary purposes
    13     only, I am including them here so that I can continue coding while our
    14     basic data types are still being implemented.
    15  *****************************************************************************/
    16 typedef enum {
    17     PS_TYPE_CHAR,
    18     PS_TYPE_SHORT,
    19     PS_TYPE_INT,
    20     PS_TYPE_LONG,
    21     PS_TYPE_UCHAR,
    22     PS_TYPE_USHORT,
    23     PS_TYPE_UINT,
    24     PS_TYPE_ULONG,
    25     PS_TYPE_FLOAT,
    26     PS_TYPE_DOUBLE,
    27     PS_TYPE_COMPLEX,
    28     PS_TYPE_OTHER,
    29 } psElemType;
    30 
    31 typedef enum {
    32     PS_DIMEN_SCALAR,
    33     PS_DIMEN_VECTOR,
    34     PS_DIMEN_TRANSV,
    35     PS_DIMEN_IMAGE,
    36     PS_DIMEN_OTHER
    37 } psDimen;
    38 
    39 typedef struct
    40 {
    41     psElemType type;
    42     psDimen dimen;
    43 }
    44 psType;
    45 
    46 typedef struct
    47 {
    48     psType type;                ///< Type of data.
    49     int nalloc;                 ///< Total number of elements available.
    50     int n;                      ///< Number of elements in use.
    51 
    52     union {
    53         int *vecI;
    54         float *vecF;
    55         double *vecD;
    56         //        complex float *vecC;
    57         void **vecP;
    58     }vec;                       ///< Union with array data.
    59 }
    60 psVector;
    61 
    62 psVector *psVectorAlloc(int        nalloc,
    63                         psElemType type)
    64 {
    65     return(NULL);
    66 }
    67 
    68 psVector *psVectorRealloc(psVector myVector,
    69                           int nalloc)
    70 {
    71     return(NULL);
    72 }
    73 
    74 void psVectorFree(psVector *restrict psArr)
    75 {}
     11int mine(psVector *myVector);
    7612
    7713/** statistics which may be calculated */
    … …  
    13369/** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
    13470psStats *
    135 psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed
    136              // must be FLOAT
    137              const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0
    138              // must be INT or NULL
     71psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed: must be FLOAT
     72             const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
    13973             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskVector
    14074             psStats *stats  ///< stats structure defines stats to be calculated and how

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

    r643 r652  
    33#include <string.h>
    44#include <stdarg.h>
    5 #include "pslib.h"
    6 /******************************************************************************
     5/******************************************************************************
     6//#include "pslib.h"
    77//#include "psMemory.h"
    88//#include "psTrace.h"
    … …  
    6464    psHistogram *newHist = NULL;
    6565    int numBins = 0;
     66    psType dummyFloatType;
     67    psType dummyIntType;
     68
     69    dummyFloatType.type = PS_TYPE_FLOAT;
     70    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     71    dummyIntType.type = PS_TYPE_FLOAT;
     72    dummyIntType.dimen = PS_DIMEN_VECTOR;
    6673
    6774    numBins = 1 + ((upper - lower) / size);
    6875    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    69     newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    70     newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    71     newHist->nums =  psVectorAlloc(numBins, PS_TYPE_INT);
     76    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
     77    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
     78    newHist->nums =  psVectorAlloc(dummyIntType, numBins);
    7279    newHist->minVal = lower;
    7380    newHist->maxVal = upper;
    … …  
    8693    int numBins = 0;
    8794    int i;
     95    psType dummyFloatType;
     96    psType dummyIntType;
     97
     98    dummyFloatType.type = PS_TYPE_FLOAT;
     99    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     100    dummyIntType.type = PS_TYPE_FLOAT;
     101    dummyIntType.dimen = PS_DIMEN_VECTOR;
    88102
    89103    if (lower->n != upper->n) {
    … …  
    92106    numBins = lower->n;
    93107    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
    94     newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
    95     newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
     108    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
     109    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
    96110    for (i=0;i<numBins;i++) {
    97         newHist->lower->vec.vecF[i] = lower->vec.vecF[i];
    98         newHist->upper->vec.vecF[i] = upper->vec.vecF[i];
    99     }
    100     newHist->nums = psVectorAlloc(numBins, PS_TYPE_INT);
     111        newHist->lower->vec.f[i] = lower->vec.f[i];
     112        newHist->upper->vec.f[i] = upper->vec.f[i];
     113    }
     114    newHist->nums = psVectorAlloc(dummyIntType, numBins);
    101115    newHist->minVal = minVal;
    102116    newHist->maxVal = maxVal;
    … …  
    135149    binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
    136150    for (i=0;i<myVector->n;i++) {
    137         if (myVector->vec.vecF[i] < myHist->minVal) {
     151        if (myVector->vec.f[i] < myHist->minVal) {
    138152            myHist->minNum++;
    139         } else if (myVector->vec.vecF[i] > myHist->maxVal) {
     153        } else if (myVector->vec.f[i] > myHist->maxVal) {
    140154            myHist->maxNum++;
    141155        } else {
    142             binNum = (int) ((myVector->vec.vecF[i] - myHist->minVal) / binSize);
    143             if ((myVector->vec.vecF[i] >= myHist->lower->vec.vecF[i]) &&
    144                     (myVector->vec.vecF[i] <= myHist->upper->vec.vecF[i])) {
    145                 myHist->nums->vec.vecI[i]++;
     156            binNum = (int) ((myVector->vec.f[i] - myHist->minVal) / binSize);
     157            if ((myVector->vec.f[i] >= myHist->lower->vec.f[i]) &&
     158                    (myVector->vec.f[i] <= myHist->upper->vec.f[i])) {
     159                myHist->nums->vec.i[i]++;
    146160            } else {
    147161                psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
    … …  
    164178    if (maskVector != NULL) {
    165179        for (i=0;i<myVector->n;i++) {
    166             if (!(maskVal & maskVector->vec.vecI[i])) {
    167                 mean+= maskVector->vec.vecI[i];
     180            if (!(maskVal & maskVector->vec.i[i])) {
     181                mean+= maskVector->vec.i[i];
    168182            }
    169183        }
    170184    } else {
    171185        for (i=0;i<myVector->n;i++) {
    172             mean+= maskVector->vec.vecI[i];
     186            mean+= maskVector->vec.i[i];
    173187        }
    174188    }
    … …  
    185199    if (maskVector != NULL) {
    186200        for (i=0;i<myVector->n;i++) {
    187             if (!(maskVal & maskVector->vec.vecI[i])) {
    188                 if (myVector->vec.vecF[i] > max) {
    189                     max = maskVector->vec.vecI[i];
     201            if (!(maskVal & maskVector->vec.i[i])) {
     202                if (myVector->vec.f[i] > max) {
     203                    max = maskVector->vec.i[i];
    190204                }
    191205            }
    … …  
    193207    } else {
    194208        for (i=0;i<myVector->n;i++) {
    195             if (myVector->vec.vecF[i] > max) {
    196                 max = maskVector->vec.vecI[i];
     209            if (myVector->vec.f[i] > max) {
     210                max = maskVector->vec.i[i];
    197211            }
    198212        }
    … …  
    210224    if (maskVector != NULL) {
    211225        for (i=0;i<myVector->n;i++) {
    212             if (!(maskVal & maskVector->vec.vecI[i])) {
    213                 if (myVector->vec.vecF[i] < min) {
    214                     min = maskVector->vec.vecI[i];
     226            if (!(maskVal & maskVector->vec.i[i])) {
     227                if (myVector->vec.f[i] < min) {
     228                    min = maskVector->vec.i[i];
    215229                }
    216230            }
    … …  
    218232    } else {
    219233        for (i=0;i<myVector->n;i++) {
    220             if (myVector->vec.vecF[i] < min) {
    221                 min = maskVector->vec.vecI[i];
     234            if (myVector->vec.f[i] < min) {
     235                min = maskVector->vec.i[i];
    222236            }
    223237        }
    … …  
    226240}
    227241
     242/******************************************************************************
     243 
     244 *****************************************************************************/
    228245int p_psArrayNValues(const psVector *restrict myVector, // float
    229246                     const psVector *restrict maskVector, // ints
    … …  
    235252    if (maskVector != NULL) {
    236253        for (i=0;i<myVector->n;i++) {
    237             if (!(maskVal & maskVector->vec.vecI[i])) {
     254            if (!(maskVal & maskVector->vec.i[i])) {
    238255                numData++;
    239256            }
    … …  
    243260    }
    244261    return(numData);
     262}
     263
     264
     265
     266#define MEDIAN_SIZE_THRESHOLD 10000
     267float p_psArraySampleMedian(const psVector *restrict myVector,
     268                            const psVector *restrict maskVector,
     269                            unsigned int maskVal)
     270{
     271    psVector *unsortedVector = NULL;
     272    psVector *sortedVector = NULL;
     273    psType dummyFloatType;
     274    int dataSize = 0;
     275    int count = 0;
     276    int i = 0;
     277    float median = 0.0;
     278
     279    dummyFloatType.type = PS_TYPE_FLOAT;
     280    dummyFloatType.dimen = PS_DIMEN_VECTOR;
     281    dataSize = p_psArrayNValues(myVector, maskVector, maskVal);
     282
     283    if (dataSize < MEDIAN_SIZE_THRESHOLD) {
     284        unsortedVector = psVectorAlloc(dummyFloatType, dataSize);
     285        sortedVector   = psVectorAlloc(dummyFloatType, dataSize);
     286
     287        count = 0;
     288        if (maskVector != NULL) {
     289            for (i=0;i<myVector->n;i++) {
     290                if (!(maskVal & maskVector->vec.i[i])) {
     291                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     292                }
     293            }
     294            psSort(sortedVector, unsortedVector);
     295        } else {
     296            psSort(sortedVector, myVector);
     297        }
     298
     299        if (0 == (dataSize % 2)) {
     300            median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
     301                            sortedVector->vec.f[dataSize/2]);
     302        } else {
     303            median = sortedVector->vec.f[dataSize/2];
     304        }
     305    } else {
     306        // BROAD: Calculate the Robust Median
     307        // Determine the LQ of the distribution.
     308        // Determine the UQ of the distribution.
     309        // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
     310        // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
     311        // Find the bin with the peak value between LQ and UQ (the MODE)
     312        // dL = (UQ - LQ) / 8
     313        // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     314        // The resulting fit parameters are the robust mean, mean_r, and sigma
     315        psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
     316    }
     317
     318    psVectorFree(unsortedVector);
     319    psVectorFree(sortedVector);
     320    return(median);
     321}
     322
     323
     324int p_psArrayXXX(const psVector *restrict myVector,
     325                 const psVector *restrict maskVector,
     326                 unsigned int maskVal)
     327{
     328    printf("ERROR: Don't call me: p_psArrayXXX()\n");
     329    exit(1);
    245330}
    246331
    … …  
    265350    switch (stats->options) {
    266351    case PS_STAT_SAMPLE_MEAN:
    267         newStruct->max = p_psArraySampleMean(myVector, maskVector, maskVal);
     352        newStruct->sampleMean = p_psArraySampleMean(myVector, maskVector, maskVal);
    268353        break;
    269354    case PS_STAT_SAMPLE_MEDIAN:
    270         printf("Code me!\n");
     355        newStruct->sampleMedian = p_psArrayXXX(myVector, maskVector, maskVal);
    271356        break;
    272357    case PS_STAT_SAMPLE_STDEV:
    273         printf("Code me!\n");
     358        newStruct->sampleStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    274359        break;
    275360    case PS_STAT_SAMPLE_UQ:
    276         printf("Code me!\n");
     361        newStruct->sampleUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    277362        break;
    278363    case PS_STAT_SAMPLE_LQ:
    279         printf("Code me!\n");
     364        newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    280365        break;
    281366    case PS_STAT_ROBUST_MEAN:
    282         printf("Code me!\n");
     367        newStruct->robustMean = p_psArrayXXX(myVector, maskVector, maskVal);
    283368        break;
    284369    case PS_STAT_ROBUST_MEAN_NVALUES:
    285         printf("Code me!\n");
     370        newStruct->robustMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    286371        break;
    287372    case PS_STAT_ROBUST_MEDIAN:
    288         printf("Code me!\n");
     373        newStruct->robustMedian = p_psArrayXXX(myVector, maskVector, maskVal);
    289374        break;
    290375    case PS_STAT_ROBUST_MEDIAN_NVALUES:
    291         printf("Code me!\n");
     376        newStruct->robustMedianNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    292377        break;
    293378    case PS_STAT_ROBUST_MODE:
    294         printf("Code me!\n");
     379        newStruct->robustMode = p_psArrayXXX(myVector, maskVector, maskVal);
    295380        break;
    296381    case PS_STAT_ROBUST_MODE_NVALUES:
    297         printf("Code me!\n");
     382        newStruct->robustModeNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    298383        break;
    299384    case PS_STAT_ROBUST_STDEV:
    300         printf("Code me!\n");
     385        newStruct->robustStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    301386        break;
    302387    case PS_STAT_ROBUST_UQ:
    303         printf("Code me!\n");
     388        newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    304389        break;
    305390    case PS_STAT_ROBUST_LQ:
    306         printf("Code me!\n");
     391        newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    307392        break;
    308393    case PS_STAT_CLIPPED_MEAN:
    309         printf("Code me!\n");
     394        newStruct->clippedMean = p_psArrayXXX(myVector, maskVector, maskVal);
    310395        break;
    311396    case PS_STAT_CLIPPED_MEAN_NVALUES:
    312         printf("Code me!\n");
     397        newStruct->clippedMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
    313398        break;
    314399    case PS_STAT_CLIPPED_MEAN_NSIGMA:
    315         printf("Code me!\n");
     400        newStruct->clipSigma = p_psArrayXXX(myVector, maskVector, maskVal);
    316401        break;
    317402    case PS_STAT_CLIPPED_STDEV:
    318         printf("Code me!\n");
     403        newStruct->clippedStdev = p_psArrayXXX(myVector, maskVector, maskVal);
    319404        break;
    320405    case PS_STAT_MAX:

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

    r644 r652  
    66 *  \ingroup MathGroup
    77 */
     8#include "psLib.h"
     9#include "psVector.h"
    810
    9 /******************************************************************************
    10     The following typedefs and functions belong in the psArray.h and psArray.c
    11     files.  However, those files are not available at this time, due to the
    12     changing definition of the psLib data/functions.  For temporary purposes
    13     only, I am including them here so that I can continue coding while our
    14     basic data types are still being implemented.
    15  *****************************************************************************/
    16 typedef enum {
    17     PS_TYPE_CHAR,
    18     PS_TYPE_SHORT,
    19     PS_TYPE_INT,
    20     PS_TYPE_LONG,
    21     PS_TYPE_UCHAR,
    22     PS_TYPE_USHORT,
    23     PS_TYPE_UINT,
    24     PS_TYPE_ULONG,
    25     PS_TYPE_FLOAT,
    26     PS_TYPE_DOUBLE,
    27     PS_TYPE_COMPLEX,
    28     PS_TYPE_OTHER,
    29 } psElemType;
    30 
    31 typedef enum {
    32     PS_DIMEN_SCALAR,
    33     PS_DIMEN_VECTOR,
    34     PS_DIMEN_TRANSV,
    35     PS_DIMEN_IMAGE,
    36     PS_DIMEN_OTHER
    37 } psDimen;
    38 
    39 typedef struct
    40 {
    41     psElemType type;
    42     psDimen dimen;
    43 }
    44 psType;
    45 
    46 typedef struct
    47 {
    48     psType type;                ///< Type of data.
    49     int nalloc;                 ///< Total number of elements available.
    50     int n;                      ///< Number of elements in use.
    51 
    52     union {
    53         int *vecI;
    54         float *vecF;
    55         double *vecD;
    56         //        complex float *vecC;
    57         void **vecP;
    58     }vec;                       ///< Union with array data.
    59 }
    60 psVector;
    61 
    62 psVector *psVectorAlloc(int        nalloc,
    63                         psElemType type)
    64 {
    65     return(NULL);
    66 }
    67 
    68 psVector *psVectorRealloc(psVector myVector,
    69                           int nalloc)
    70 {
    71     return(NULL);
    72 }
    73 
    74 void psVectorFree(psVector *restrict psArr)
    75 {}
     11int mine(psVector *myVector);
    7612
    7713/** statistics which may be calculated */
    … …  
    13369/** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
    13470psStats *
    135 psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed
    136              // must be FLOAT
    137              const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0
    138              // must be INT or NULL
     71psArrayStats(const psVector *restrict myVector, ///< Vector to be analysed: must be FLOAT
     72             const psVector *restrict maskVector, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
    13973             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskVector
    14074             psStats *stats  ///< stats structure defines stats to be calculated and how
Note: See TracChangeset for help on using the changeset viewer.