IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 26, 2009, 8:40:07 PM (17 years ago)
Author:
eugene
Message:

incorporating changes from 16bit mask upgrades (eam_branch_20081230)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/math/psStats.c

    r20515 r21183  
    1313 * use ->min and ->max (PS_STAT_USE_RANGE)
    1414 *
    15  *  @version $Revision: 1.230 $ $Name: not supported by cvs2svn $
    16  *  @date $Date: 2008-11-04 00:55:14 $
     15 *  @version $Revision: 1.231 $ $Name: not supported by cvs2svn $
     16 *  @date $Date: 2009-01-27 06:39:38 $
    1717 *
    1818 *  Copyright 2006 IfA, University of Hawaii
     
    6868#define PS_CLIPPED_SIGMA_UB 10.0
    6969#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
    70 #define MASK_MARK 0x80   // bit to use internally to mark data as bad
     70
     71#define MASK_MARK 0x80   // XXX : can we change this? bit to use internally to mark data as bad
    7172#define PS_ROBUST_MAX_ITERATIONS 20     // Maximum number of iterations for robust statistics
    7273
     
    178179                                 const psVector* errors,
    179180                                 const psVector* maskVector,
    180                                  psMaskType maskVal,
     181                                 psVectorMaskType maskVal,
    181182                                 psStats* stats)
    182183{
     
    190191    int numData = myVector->n;          // Number of data points
    191192
    192     psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
     193    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
    193194    bool useRange = stats->options & PS_STAT_USE_RANGE;
    194195
     
    245246    static long vectorMinMax(const psVector* myVector,
    246247                             const psVector* maskVector,
    247                              psMaskType maskVal,
     248                             psVectorMaskType maskVal,
    248249                             psStats* stats
    249250        )
     
    257258    int numValid = 0;                   // Number of valid values
    258259
    259     psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
     260    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
    260261    bool useRange = stats->options & PS_STAT_USE_RANGE;
    261262
     
    302303static bool vectorSampleMedian(const psVector* inVector,
    303304                               const psVector* maskVector,
    304                                psMaskType maskVal,
     305                               psVectorMaskType maskVal,
    305306                               psStats* stats)
    306307{
     
    308309
    309310    bool useRange = stats->options & PS_STAT_USE_RANGE;
    310     psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector
     311    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector
    311312    psF32 *input = inVector->data.F32; // Dereference the vector
    312313
     
    387388                              const psVector* errors,
    388389                              const psVector* maskVector,
    389                               psMaskType maskVal,
     390                              psVectorMaskType maskVal,
    390391                              psStats* stats)
    391392{
     
    407408
    408409    psF32 *data = myVector->data.F32;   // Dereference
    409     psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
     410    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
    410411    bool useRange = stats->options & PS_STAT_USE_RANGE;
    411412    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;
     
    468469static bool vectorSampleMoments(const psVector* myVector,
    469470                                const psVector* maskVector,
    470                                 psMaskType maskVal,
     471                                psVectorMaskType maskVal,
    471472                                psStats* stats)
    472473{
     
    490491
    491492    psF32 *data = myVector->data.F32;   // Dereference
    492     psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
     493    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA;
    493494    bool useRange = stats->options & PS_STAT_USE_RANGE;
    494495
     
    557558                               const psVector* errors,
    558559                               psVector* maskInput,
    559                                psMaskType maskValInput,
     560                               psVectorMaskType maskValInput,
    560561                               psStats* stats
    561562    )
     
    580581
    581582    // We copy the mask vector, to preserve the original
    582     psMaskType maskVal = MASK_MARK | maskValInput;
     583    psVectorMaskType maskVal = MASK_MARK | maskValInput;
    583584
    584585    // use the temporary vector for local temporary mask
    585     stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_U8);
     586    stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK);
    586587    psVector *tmpMask = stats->tmpMask;
    587588    psVectorInit(tmpMask, 0);
    588589    if (maskInput) {
    589590        for (long i = 0; i < myVector->n; i++) {
    590             if (maskInput->data.U8[i] & maskValInput) {
    591                 tmpMask->data.U8[i] = maskVal;
     591            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) {
     592                tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
    592593            }
    593594        }
     
    629630            // sqrt(A))
    630631            for (long j = 0; j < myVector->n; j++) {
    631                 if (!tmpMask->data.U8[j] &&
     632                if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] &&
    632633                    fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
    633                     tmpMask->data.U8[j] = 0xff;
     634                    tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
    634635                    psTrace(TRACE, 10, "Clipped %ld: %f +/- %f\n", j,
    635636                            myVector->data.F32[j], errors->data.F32[j]);
     
    640641        } else {
    641642            for (long j = 0; j < myVector->n; j++) {
    642                 if (!tmpMask->data.U8[j] &&
     643                if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] &&
    643644                    fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
    644                     tmpMask->data.U8[j] = 0xff;
     645                    tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff;
    645646                    psTrace(TRACE, 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]);
    646647                    numClipped++;
     
    715716                              const psVector* errors,
    716717                              psVector* maskInput,
    717                               psMaskType maskValInput,
     718                              psVectorMaskType maskValInput,
    718719                              psStats* stats)
    719720{
     
    726727    // and tested even if there is no supplied mask (and/or the maskVal is 0)
    727728    // XXX this would be better if we had globally defined mask values
    728     psMaskType maskVal = MASK_MARK | maskValInput;
    729     psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_MASK); // The actual mask we will use
     729    psVectorMaskType maskVal = MASK_MARK | maskValInput;
     730    psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_VECTOR_MASK); // The actual mask we will use
    730731    psVectorInit(mask, 0);
    731732    if (maskInput) {
    732733        for (long i = 0; i < myVector->n; i++) {
    733             if (maskInput->data.U8[i] & maskValInput) {
    734                 mask->data.U8[i] = maskVal;
     734            if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) {
     735                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
    735736            }
    736737        }
     
    934935                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
    935936                    // XXXX is this correct?  is MASK_MARK safe?
    936                     mask->data.U8[i] |= MASK_MARK;
     937                    mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
    937938                    psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
    938939                }
     
    998999    long N50 = 0;
    9991000    for (long i = 0 ; i < myVector->n ; i++) {
    1000         if (!mask->data.U8[i] &&
     1001        if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[i] &&
    10011002            (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
    10021003            N50++;
     
    10461047                               const psVector* errors,
    10471048                               psVector* mask,
    1048                                psMaskType maskVal,
     1049                               psVectorMaskType maskVal,
    10491050                               psStats* stats)
    10501051{
     
    12241225                                  const psVector* errors,
    12251226                                  psVector* mask,
    1226                                   psMaskType maskVal,
     1227                                  psVectorMaskType maskVal,
    12271228                                  psStats* stats)
    12281229{
     
    13471348        // fitStats->clipIter = 3.0;
    13481349        // fitStats->clipSigma = 3.0;
    1349         // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_U8);
     1350        // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_VECTOR_MASK);
    13501351        // psVectorInit (fitMask, 0);
    13511352
     
    14181419                                  const psVector* errors,
    14191420                                  psVector* mask,
    1420                                   psMaskType maskVal,
     1421                                  psVectorMaskType maskVal,
    14211422                                  psStats* stats)
    14221423{
     
    17131714                                  const psVector* errors,
    17141715                                  psVector* mask,
    1715                                   psMaskType maskVal,
     1716                                  psVectorMaskType maskVal,
    17161717                                  psStats* stats)
    17171718{
     
    22252226                   const psVector* errors,
    22262227                   const psVector* mask,
    2227                    psMaskType maskVal)
     2228                   psVectorMaskType maskVal)
    22282229{
    22292230    psTrace(TRACE, 3,"---- %s() begin  ----\n", __func__);
     
    22332234    if (mask) {
    22342235        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, false);
    2235         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false);
     2236        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
    22362237    }
    22372238    if (errors) {
     
    22552256        }
    22562257    }
    2257     psVector *maskU8 = NULL;            // Input mask vector, U8 version
     2258    psVector *maskVector = NULL;            // Input mask vector, U8 version
    22582259    if (mask) {
    2259         if (mask->type.type == PS_TYPE_MASK) {
    2260             maskU8 = psMemIncrRefCounter((psPtr)mask);
     2260        if (mask->type.type == PS_TYPE_VECTOR_MASK) {
     2261            maskVector = psMemIncrRefCounter((psPtr)mask);
    22612262        } else {
    2262             maskU8 = psVectorCopy(NULL, mask, PS_TYPE_MASK);
     2263            maskVector = psVectorCopy(NULL, mask, PS_TYPE_VECTOR_MASK);
    22632264        }
    22642265    }
     
    22792280    // ************************************************************************
    22802281    if (stats->options & PS_STAT_SAMPLE_MEAN) {
    2281         if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
     2282        if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats)) {
    22822283            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean");
    22832284            status &= false;
     
    22872288    // ************************************************************************
    22882289    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
    2289         if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
     2290        if (!vectorSampleMedian(inF32, maskVector, maskVal, stats)) {
    22902291            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
    22912292            status &= false;
     
    22952296    // ************************************************************************
    22962297    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    2297         if (!vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats)) {
     2298        if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats)) {
    22982299            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev");
    22992300            status &= false;
     
    23022303
    23032304    if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) {
    2304         if (!vectorSampleMoments(inF32, maskU8, maskVal, stats)) {
     2305        if (!vectorSampleMoments(inF32, maskVector, maskVal, stats)) {
    23052306            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample moments");
    23062307            status &= false;
     
    23102311    // ************************************************************************
    23112312    if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) {
    2312         if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) {
     2313        if (vectorMinMax(inF32, maskVector, maskVal, stats) == 0) {
    23132314            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum");
    23142315            status &= false;
     
    23182319    // ************************************************************************
    23192320    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) {
    2320         if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
     2321        if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats)) {
    23212322            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics"));
    23222323            status &= false;
     
    23262327    // ************************************************************************
    23272328    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
    2328         if (!vectorFittedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
     2329        if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
    23292330            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    23302331            status &= false;
     
    23372338            psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V2");
    23382339        }
    2339         if (!vectorFittedStats_v2(inF32, errorsF32, maskU8, maskVal, stats)) {
     2340        if (!vectorFittedStats_v2(inF32, errorsF32, maskVector, maskVal, stats)) {
    23402341            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    23412342            status &= false;
     
    23482349            psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V3");
    23492350        }
    2350         if (!vectorFittedStats_v3(inF32, errorsF32, maskU8, maskVal, stats)) {
     2351        if (!vectorFittedStats_v3(inF32, errorsF32, maskVector, maskVal, stats)) {
    23512352            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    23522353            status &= false;
     
    23592360            psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4");
    23602361        }
    2361         if (!vectorFittedStats_v4(inF32, errorsF32, maskU8, maskVal, stats)) {
     2362        if (!vectorFittedStats_v4(inF32, errorsF32, maskVector, maskVal, stats)) {
    23622363            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    23632364            status &= false;
     
    23672368    // ************************************************************************
    23682369    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    2369         if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
     2370        if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
    23702371            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n");
    23712372            status &= false;
     
    23752376    psFree(inF32);
    23762377    psFree(errorsF32);
    2377     psFree(maskU8);
     2378    psFree(maskVector);
    23782379    psTrace(TRACE, 3,"---- %s() end  ----\n", __func__);
    23792380    return status;
Note: See TracChangeset for help on using the changeset viewer.