IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 699


Ignore:
Timestamp:
May 14, 2004, 4:02:22 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r683 r699  
    270270    psVector *unsortedVector = NULL;
    271271    psVector *sortedVector = NULL;
    272     int dataSize = 0;
    273272    int count = 0;
    274273    int i = 0;
    275274    float median = 0.0;
    276275
    277     p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
    278 
    279     if (dataSize < MEDIAN_SIZE_THRESHOLD) {
    280         unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
    281         sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
    282 
    283         count = 0;
    284         if (maskVector != NULL) {
    285             for (i=0;i<myVector->n;i++) {
    286                 if (!(maskVal & maskVector->vec.ui8[i])) {
    287                     unsortedVector->vec.f[count++] = maskVector->vec.f[i];
    288                 }
    289             }
    290             psSort(sortedVector, unsortedVector);
    291         } else {
    292             psSort(sortedVector, myVector);
    293         }
    294 
    295         if (0 == (dataSize % 2)) {
    296             median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
    297                             sortedVector->vec.f[dataSize/2]);
    298         } else {
    299             median = sortedVector->vec.f[dataSize/2];
    300         }
    301     } else {
    302         // BROAD: Calculate the Robust Median
    303         // Determine the LQ of the distribution.
    304         // Determine the UQ of the distribution.
    305         // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
    306         // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
    307         // Find the bin with the peak value between LQ and UQ (the MODE)
    308         // dL = (UQ - LQ) / 8
    309         // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    310         // The resulting fit parameters are the robust mean, mean_r, and sigma
    311         psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
     276    if (-1 == newStruct->nValues) {
     277        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     278    }
     279    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     280    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     281
     282    if (maskVector != NULL) {
     283        for (i=0;i<myVector->n;i++) {
     284            if (!(maskVal & maskVector->vec.ui8[i])) {
     285                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     286            }
     287        }
     288        psSort(sortedVector, unsortedVector);
     289    } else {
     290        psSort(sortedVector, myVector);
     291    }
     292
     293    if (0 == (newStruct->nValues % 2)) {
     294        median = 0.5 * (sortedVector->vec.f[(newStruct->nValues/2)-1] +
     295                        sortedVector->vec.f[newStruct->nValues/2]);
     296    } else {
     297        median = sortedVector->vec.f[newStruct->nValues/2];
    312298    }
    313299
    … …  
    315301    psVectorFree(sortedVector);
    316302    newStruct->sampleMedian = median;
     303}
     304
     305void p_psArrayRobustMedian(const psVector *restrict myVector,
     306                           const psVector *restrict maskVector,
     307                           unsigned int maskVal,
     308                           psStats *newStruct)
     309{
     310    psHistogram *robustHistogram = NULL;
     311    float binSize = 0.0;
     312
     313    //    if (isnan(myVector->robustLQ) ||
     314    //        isnan(myVector->robustUQ)) {
     315    //        p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct);
     316    //    }
     317    //    binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0;
     318
     319    robustHistogram = psHistogramAlloc(newStruct->min,
     320                                       newStruct->max,
     321                                       binSize);
     322    //    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     323    //    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
     324    //    p_psArraySmooth(robustHistogram, (binSize / 4.0));
     325    //    dL = (myVector->robustUQ - myVector->robustLQ) / 8.0;
     326
     327
     328    // BROAD: Calculate the Robust Median
     329    // Determine the LQ of the distribution.
     330    // Determine the UQ of the distribution.
     331    // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
     332    // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
     333
     334    // Find the bin with the peak value between LQ and UQ (the MODE)
     335    // dL = (UQ - LQ) / 8
     336    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     337    // The resulting fit parameters are the robust mean, mean_r, and sigma
    317338}
    318339
    … …  
    368389}
    369390
    370 void p_psArraySampleUQ(const psVector *restrict myVector,
    371                        const psVector *restrict maskVector,
    372                        unsigned int maskVal,
    373                        psStats *newStruct)
    374 {
    375     psHistogram    *nonRobustHistogram = NULL;
     391/******************************************************************************
     392 *****************************************************************************/
     393void p_psArraySampleQuartiles(const psVector *restrict myVector,
     394                              const psVector *restrict maskVector,
     395                              unsigned int maskVal,
     396                              psStats *newStruct)
     397{
     398    psVector *unsortedVector = NULL;
     399    psVector *sortedVector = NULL;
     400    int count = 0;
     401    int ind = 0;
     402    int i = 0;
     403
     404    // return is we have already calculated both quartile points.
     405    if ((!isnan(newStruct->sampleLQ)) &&
     406            (!isnan(newStruct->sampleUQ))) {
     407        return;
     408    }
    376409
    377410    if (-1 == newStruct->nValues) {
    … …  
    379412    }
    380413
    381     if (newStruct->nValues < MEDIAN_SIZE_THRESHOLD) {
    382 
    383     }
    384     else {
    385         if (0 != isnan(newStruct->sampleStdev)) {
    386             p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
    387         }
    388         if (0 != isnan(newStruct->min)) {
    389             p_psArrayMin(myVector, maskVector, maskVal, newStruct);
    390         }
    391         if (0 != isnan(newStruct->max)) {
    392             p_psArrayMax(myVector, maskVector, maskVal, newStruct);
    393         }
    394 
    395         nonRobustHistogram = psHistogramAlloc(newStruct->min,
    396                                               newStruct->max,
    397                                               10);
    398     }
    399 }
    400 
     414
     415    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     416    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     417
     418    count = 0;
     419    if (maskVector != NULL) {
     420        for (i=0;i<myVector->n;i++) {
     421            if (!(maskVal & maskVector->vec.ui8[i])) {
     422                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     423            }
     424        }
     425        psSort(sortedVector, unsortedVector);
     426    } else {
     427        psSort(sortedVector, myVector);
     428    }
     429
     430    ind = 3 * (newStruct->nValues / 4);
     431    newStruct->sampleUQ = sortedVector->vec.f[ind];
     432    ind = (newStruct->nValues / 4);
     433    newStruct->sampleLQ = sortedVector->vec.f[ind];
     434
     435    psVectorFree(unsortedVector);
     436    psVectorFree(sortedVector);
     437}
    401438
    402439/******************************************************************************
    … …  
    436473    newStruct = psStatsAlloc(stats->options);
    437474
     475    // ************************************************************************
    438476    if (stats->options & PS_STAT_SAMPLE_MEAN) {
    439477        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
    440478    }
    441479
     480    // ************************************************************************
    442481    if (stats->options & PS_STAT_MAX) {
    443482        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
    444483    }
    445484
     485    // ************************************************************************
    446486    if (stats->options & PS_STAT_MIN) {
    447487        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
    448488    }
    449489
     490    // ************************************************************************
    450491    if (stats->options & PS_STAT_NVALUES) {
    451492        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
    452493    }
    453494
    454     if (stats->options & PS_STAT_SAMPLE_MEDIAN) {}
    455 
     495    // ************************************************************************
     496    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
     497        p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
     498    }
     499
     500    // ************************************************************************
    456501    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    457502        p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
    458503    }
    459504
    460     if (stats->options & PS_STAT_SAMPLE_UQ) {
    461         p_psArraySampleUQ(myVector, maskVector, maskVal, newStruct);
    462     }
    463 
    464     if (stats->options & PS_STAT_SAMPLE_LQ) {
    465         newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    466     }
     505    // ************************************************************************
     506    if ((stats->options & PS_STAT_SAMPLE_UQ) ||
     507            (stats->options & PS_STAT_SAMPLE_LQ)) {
     508        p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct);
     509    }
     510
     511
     512
    467513
    468514    if (stats->options & PS_STAT_ROBUST_MEAN) {
    … …  
    494540    }
    495541
    496     if (stats->options & PS_STAT_ROBUST_UQ) {
    497         newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    498     }
    499 
    500     if (stats->options & PS_STAT_ROBUST_LQ) {
     542    if ((stats->options & PS_STAT_ROBUST_UQ) ||
     543            (stats->options & PS_STAT_ROBUST_LQ)) {
    501544        newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    502545    }

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

    r683 r699  
    270270    psVector *unsortedVector = NULL;
    271271    psVector *sortedVector = NULL;
    272     int dataSize = 0;
    273272    int count = 0;
    274273    int i = 0;
    275274    float median = 0.0;
    276275
    277     p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
    278 
    279     if (dataSize < MEDIAN_SIZE_THRESHOLD) {
    280         unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
    281         sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, dataSize);
    282 
    283         count = 0;
    284         if (maskVector != NULL) {
    285             for (i=0;i<myVector->n;i++) {
    286                 if (!(maskVal & maskVector->vec.ui8[i])) {
    287                     unsortedVector->vec.f[count++] = maskVector->vec.f[i];
    288                 }
    289             }
    290             psSort(sortedVector, unsortedVector);
    291         } else {
    292             psSort(sortedVector, myVector);
    293         }
    294 
    295         if (0 == (dataSize % 2)) {
    296             median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
    297                             sortedVector->vec.f[dataSize/2]);
    298         } else {
    299             median = sortedVector->vec.f[dataSize/2];
    300         }
    301     } else {
    302         // BROAD: Calculate the Robust Median
    303         // Determine the LQ of the distribution.
    304         // Determine the UQ of the distribution.
    305         // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
    306         // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
    307         // Find the bin with the peak value between LQ and UQ (the MODE)
    308         // dL = (UQ - LQ) / 8
    309         // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    310         // The resulting fit parameters are the robust mean, mean_r, and sigma
    311         psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
     276    if (-1 == newStruct->nValues) {
     277        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     278    }
     279    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     280    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     281
     282    if (maskVector != NULL) {
     283        for (i=0;i<myVector->n;i++) {
     284            if (!(maskVal & maskVector->vec.ui8[i])) {
     285                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     286            }
     287        }
     288        psSort(sortedVector, unsortedVector);
     289    } else {
     290        psSort(sortedVector, myVector);
     291    }
     292
     293    if (0 == (newStruct->nValues % 2)) {
     294        median = 0.5 * (sortedVector->vec.f[(newStruct->nValues/2)-1] +
     295                        sortedVector->vec.f[newStruct->nValues/2]);
     296    } else {
     297        median = sortedVector->vec.f[newStruct->nValues/2];
    312298    }
    313299
    … …  
    315301    psVectorFree(sortedVector);
    316302    newStruct->sampleMedian = median;
     303}
     304
     305void p_psArrayRobustMedian(const psVector *restrict myVector,
     306                           const psVector *restrict maskVector,
     307                           unsigned int maskVal,
     308                           psStats *newStruct)
     309{
     310    psHistogram *robustHistogram = NULL;
     311    float binSize = 0.0;
     312
     313    //    if (isnan(myVector->robustLQ) ||
     314    //        isnan(myVector->robustUQ)) {
     315    //        p_psArrayRobustQuartiles(myVector, maskVector, maskVal, newStruct);
     316    //    }
     317    //    binSize = ((myVector->robustUQ - myVector->robustLQ) / 1.34) / 10.0;
     318
     319    robustHistogram = psHistogramAlloc(newStruct->min,
     320                                       newStruct->max,
     321                                       binSize);
     322    //    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
     323    //    robustHistogram = psGetArrayHistogram(robustHistogram, myVector);
     324    //    p_psArraySmooth(robustHistogram, (binSize / 4.0));
     325    //    dL = (myVector->robustUQ - myVector->robustLQ) / 8.0;
     326
     327
     328    // BROAD: Calculate the Robust Median
     329    // Determine the LQ of the distribution.
     330    // Determine the UQ of the distribution.
     331    // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
     332    // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
     333
     334    // Find the bin with the peak value between LQ and UQ (the MODE)
     335    // dL = (UQ - LQ) / 8
     336    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
     337    // The resulting fit parameters are the robust mean, mean_r, and sigma
    317338}
    318339
    … …  
    368389}
    369390
    370 void p_psArraySampleUQ(const psVector *restrict myVector,
    371                        const psVector *restrict maskVector,
    372                        unsigned int maskVal,
    373                        psStats *newStruct)
    374 {
    375     psHistogram    *nonRobustHistogram = NULL;
     391/******************************************************************************
     392 *****************************************************************************/
     393void p_psArraySampleQuartiles(const psVector *restrict myVector,
     394                              const psVector *restrict maskVector,
     395                              unsigned int maskVal,
     396                              psStats *newStruct)
     397{
     398    psVector *unsortedVector = NULL;
     399    psVector *sortedVector = NULL;
     400    int count = 0;
     401    int ind = 0;
     402    int i = 0;
     403
     404    // return is we have already calculated both quartile points.
     405    if ((!isnan(newStruct->sampleLQ)) &&
     406            (!isnan(newStruct->sampleUQ))) {
     407        return;
     408    }
    376409
    377410    if (-1 == newStruct->nValues) {
    … …  
    379412    }
    380413
    381     if (newStruct->nValues < MEDIAN_SIZE_THRESHOLD) {
    382 
    383     }
    384     else {
    385         if (0 != isnan(newStruct->sampleStdev)) {
    386             p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
    387         }
    388         if (0 != isnan(newStruct->min)) {
    389             p_psArrayMin(myVector, maskVector, maskVal, newStruct);
    390         }
    391         if (0 != isnan(newStruct->max)) {
    392             p_psArrayMax(myVector, maskVector, maskVal, newStruct);
    393         }
    394 
    395         nonRobustHistogram = psHistogramAlloc(newStruct->min,
    396                                               newStruct->max,
    397                                               10);
    398     }
    399 }
    400 
     414
     415    unsortedVector = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     416    sortedVector   = psVectorAlloc(PS_TYPE_FLOAT, newStruct->nValues);
     417
     418    count = 0;
     419    if (maskVector != NULL) {
     420        for (i=0;i<myVector->n;i++) {
     421            if (!(maskVal & maskVector->vec.ui8[i])) {
     422                unsortedVector->vec.f[count++] = maskVector->vec.f[i];
     423            }
     424        }
     425        psSort(sortedVector, unsortedVector);
     426    } else {
     427        psSort(sortedVector, myVector);
     428    }
     429
     430    ind = 3 * (newStruct->nValues / 4);
     431    newStruct->sampleUQ = sortedVector->vec.f[ind];
     432    ind = (newStruct->nValues / 4);
     433    newStruct->sampleLQ = sortedVector->vec.f[ind];
     434
     435    psVectorFree(unsortedVector);
     436    psVectorFree(sortedVector);
     437}
    401438
    402439/******************************************************************************
    … …  
    436473    newStruct = psStatsAlloc(stats->options);
    437474
     475    // ************************************************************************
    438476    if (stats->options & PS_STAT_SAMPLE_MEAN) {
    439477        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
    440478    }
    441479
     480    // ************************************************************************
    442481    if (stats->options & PS_STAT_MAX) {
    443482        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
    444483    }
    445484
     485    // ************************************************************************
    446486    if (stats->options & PS_STAT_MIN) {
    447487        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
    448488    }
    449489
     490    // ************************************************************************
    450491    if (stats->options & PS_STAT_NVALUES) {
    451492        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
    452493    }
    453494
    454     if (stats->options & PS_STAT_SAMPLE_MEDIAN) {}
    455 
     495    // ************************************************************************
     496    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
     497        p_psArraySampleMedian(myVector, maskVector, maskVal, newStruct);
     498    }
     499
     500    // ************************************************************************
    456501    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    457502        p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
    458503    }
    459504
    460     if (stats->options & PS_STAT_SAMPLE_UQ) {
    461         p_psArraySampleUQ(myVector, maskVector, maskVal, newStruct);
    462     }
    463 
    464     if (stats->options & PS_STAT_SAMPLE_LQ) {
    465         newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    466     }
     505    // ************************************************************************
     506    if ((stats->options & PS_STAT_SAMPLE_UQ) ||
     507            (stats->options & PS_STAT_SAMPLE_LQ)) {
     508        p_psArraySampleQuartiles(myVector, maskVector, maskVal, newStruct);
     509    }
     510
     511
     512
    467513
    468514    if (stats->options & PS_STAT_ROBUST_MEAN) {
    … …  
    494540    }
    495541
    496     if (stats->options & PS_STAT_ROBUST_UQ) {
    497         newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
    498     }
    499 
    500     if (stats->options & PS_STAT_ROBUST_LQ) {
     542    if ((stats->options & PS_STAT_ROBUST_UQ) ||
     543            (stats->options & PS_STAT_ROBUST_LQ)) {
    501544        newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
    502545    }
Note: See TracChangeset for help on using the changeset viewer.