IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Added the median calculation functions.

File:
1 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:
Note: See TracChangeset for help on using the changeset viewer.