IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 13, 2007, 4:43:47 PM (19 years ago)
Author:
magnier
Message:

added macros PS_BIN_FOR_VALUE and PS_BIN_INTERPOLATE to clean up repeated code
converted p_psVectorBinDisect to new API
in robust stats:

  • added escape clause for consistent exit error handling
  • set PS_STATS_ROBUST_MEDIAN flag on error condition
  • forced interpolation of histogram bin value to land within bin
  • using the +/- 0.5 sigma positions rather than +/- 1.0

in fitted stats (all versions):

  • force the fitting range (binMin - binMax) to land within histogram
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r11760 r12437  
    1313 * use ->min and ->max (PS_STAT_USE_RANGE)
    1414 *
    15  *  @version $Revision: 1.201 $ $Name: not supported by cvs2svn $
    16  *  @date $Date: 2007-02-13 03:08:17 $
     15 *  @version $Revision: 1.202 $ $Name: not supported by cvs2svn $
     16 *  @date $Date: 2007-03-14 02:43:47 $
    1717 *
    1818 *  Copyright 2006 IfA, University of Hawaii
    … …  
    6767(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
    6868
     69// set the bin closest to the corresponding value.  if USE_END is +/- 1,
     70// out-of-range saturates on lower/upper bin REGARDLESS of actual value
     71#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) \
     72        tmpScalar.data.F32 = (VALUE); \
     73        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
     74        switch (result) { \
     75          case PS_BINARY_DISECT_PASS: \
     76            break; \
     77          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
     78            psTrace ("psLib.math", 6, "selected bin outside range"); \
     79            if (USE_END == -1) { RESULT = 0; } \
     80            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
     81            break; \
     82          case PS_BINARY_DISECT_INVALID_INPUT: \
     83          case PS_BINARY_DISECT_INVALID_TYPE: \
     84            psAbort ("programming error"); \
     85            break; \
     86        }
     87
     88# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
     89        float dX, dY, Xo, Yo, Xt; \
     90        if (BIN == BOUNDS->n - 1) { \
     91            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
     92            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \
     93            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     94            Yo = VECTOR->data.F32[BIN]; \
     95        } else { \
     96            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
     97            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \
     98            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     99            Yo = VECTOR->data.F32[BIN]; \
     100        } \
     101        if (dY != 0) { \
     102            Xt = (VALUE - Yo)*dX/dY + Xo; \
     103        } else { \
     104            Xt = Xo; \
     105        } \
     106        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
     107        psTrace("psLib.math", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
     108                Xo, Yo, dX, dY, Xt, VALUE); \
     109        RESULT = Xt; }
     110
    69111/*****************************************************************************/
    70112/* TYPE DEFINITIONS                                                          */
    … …  
    120162To optmize this, use a macro and ifdef in or out the three states (errors, mask, range)
    121163*****************************************************************************/
    122 static bool vectorSampleMean(const psVector* myVector,
    123                             const psVector* errors,
    124                             const psVector* maskVector,
    125                             psMaskType maskVal,
    126                             psStats* stats)
     164    static bool vectorSampleMean(const psVector* myVector,
     165                                const psVector* errors,
     166                                const psVector* maskVector,
     167                                psMaskType maskVal,
     168                                psStats* stats)
    127169{
    128170    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    188230(mask: 1, range: 0):  0.200 sec  0.244 sec
    189231*****************************************************************************/
    190 static long vectorMinMax(const psVector* myVector,
    191                          const psVector* maskVector,
    192                          psMaskType maskVal,
    193                          psStats* stats
    194                         )
     232    static long vectorMinMax(const psVector* myVector,
     233                             const psVector* maskVector,
     234                             psMaskType maskVal,
     235                             psStats* stats
     236        )
    195237{
    196238    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    331373
    332374*****************************************************************************/
    333 static bool vectorSampleStdev(const psVector* myVector,
    334                               const psVector* errors,
    335                               const psVector* maskVector,
    336                               psMaskType maskVal,
    337                               psStats* stats)
     375    static bool vectorSampleStdev(const psVector* myVector,
     376                                  const psVector* errors,
     377                                  const psVector* maskVector,
     378                                  psMaskType maskVal,
     379                                  psStats* stats)
    338380{
    339381    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    422464                               psMaskType maskValInput,
    423465                               psStats* stats
    424                               )
     466    )
    425467{
    426468    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    491533            for (long j = 0; j < myVector->n; j++) {
    492534                if (!tmpMask->data.U8[j] &&
    493                         fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
     535                    fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
    494536                    tmpMask->data.U8[j] = 0xff;
    495537                    psTrace("psLib.math", 10, "Clipped %ld: %f +/- %f\n", j,
    … …  
    502544            for (long j = 0; j < myVector->n; j++) {
    503545                if (!tmpMask->data.U8[j] &&
    504                         fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
     546                    fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
    505547                    tmpMask->data.U8[j] = 0xff;
    506548                    psTrace("psLib.math", 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]);
    … …  
    584626    }
    585627
    586     psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
     628    psVectorBinaryDisectResult result;
     629    psScalar tmpScalar;                 // Temporary scalar variable, for psVectorBinaryDisect
    587630    tmpScalar.type.type = PS_TYPE_F32;
    588631
    … …  
    609652    long totalDataPoints = 0;           // Total number of (unmasked) data points
    610653
     654    float binSize = 0.0;            // Size of bins for histogram
     655    long binLo, binHi;
     656    long binL2, binH2;
     657    long binMedian;
     658
    611659    // Iterate to get the best bin size
    612660    for (int iterate = 1; iterate > 0; iterate++) {
    … …  
    620668            // Data range calculation failed
    621669            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
    622             psFree(statsMinMax);
    623             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    624             psFree(mask);
    625             return false;
     670            goto escape;
    626671        }
    627672        psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
    … …  
    629674        // If all data points have the same value, then we set the appropriate members of stats and return.
    630675        if (fabs(max - min) <= FLT_EPSILON) {
    631             psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
    632676            stats->robustMedian = min;
    633677            stats->robustStdev = 0.0;
    … …  
    635679            stats->robustLQ = min;
    636680            stats->robustN50 = numValid;
    637             psFree(statsMinMax);
    638 
    639681            stats->results |= PS_STAT_ROBUST_MEDIAN;
    640682            stats->results |= PS_STAT_ROBUST_STDEV;
    641683            stats->results |= PS_STAT_ROBUST_QUARTILE;
    642 
     684            psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
    643685            psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
    644686            psFree(mask);
     687            psFree(statsMinMax);
    645688            return true;
    646689        }
    647690
    648         float binSize = 0.0;            // Size of bins for histogram
    649691        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
    650692            // Set initial bin size to the specified value.
    … …  
    668710        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    669711            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
    670             psFree(histogram);
    671             psFree(statsMinMax);
    672             return false;
     712            goto escape;
    673713        }
    674714        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    691731        totalDataPoints = cumulative->nums->data.F32[numBins - 1];
    692732        psTrace("psLib.math", 6, "Total data points is %ld\n", totalDataPoints);
    693         long binMedian;
    694         if (totalDataPoints/2.0 < cumulative->nums->data.F32[0]) {
    695             // Special case: the median is in the first bin.  Perhaps we should recode this.
    696             // XXX: Must try a special case where the median in in the last bin.
    697             binMedian = 0;
    698         } else {
    699             tmpScalar.data.F32 = totalDataPoints/2.0;
    700             binMedian = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    701             if (binMedian < 0) {
    702                 psError(PS_ERR_UNKNOWN, false,
    703                         "Failed to calculate the 50 precent data point (%ld).\n", binMedian);
    704                 psFree(statsMinMax);
    705                 psFree(histogram);
    706                 psFree(cumulative);
    707                 psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    708                 psFree(mask);
    709                 return false;
    710             }
    711         }
     733
     734        PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
    712735        psTrace("psLib.math", 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian,
    713736                cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
    … …  
    715738        // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
    716739        stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
    717                               binMedian, totalDataPoints/2.0);
     740                                                                binMedian, totalDataPoints/2.0);
    718741        if (isnan(stats->robustMedian)) {
    719             psError(PS_ERR_UNKNOWN, false,
    720                     "Failed to fit a quadratic and calculate the 50-percent position.\n");
    721             psFree(statsMinMax);
    722             psFree(histogram);
    723             psFree(cumulative);
    724             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    725             psFree(mask);
    726             return false;
     742            psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n");
     743            goto escape;
    727744        }
    728745        psTrace("psLib.math", 6, "Current robust median is %f\n", stats->robustMedian);
    729746
    730         // ADD step 4: Find the bins which contains the 15.8655% and 84.1345% data points.
    731         long binLo = 0;                 // Bin containing the 15.8655% mark
    732         tmpScalar.data.F32 = totalDataPoints * 0.158655f;
    733         if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
    734             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
    735             binLo = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    736             if (binLo > cumulative->nums->n-1) {
    737                 binLo = cumulative->nums->n-1;
    738             }
    739         }
    740 
    741         long binHi = 0;                 // Bin containing the 84.1345% mark
    742         tmpScalar.data.F32 = totalDataPoints * 0.841345f;
    743         if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
    744             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
    745             binHi = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    746             if (binHi > cumulative->nums->n-1) {
    747                 binHi = cumulative->nums->n-1;
    748             }
    749         }
     747        // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data points
     748        // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points
     749        PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1);
     750        PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1);
     751
     752        PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1);
     753        PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1);
     754
    750755        psTrace("psLib.math", 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
    751756                binLo, binHi);
    … …  
    755760        if ((binLo < 0) || (binHi < 0)) {
    756761            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
    757             psFree(statsMinMax);
    758             psFree(histogram);
    759             psFree(cumulative);
    760             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    761             psFree(mask);
    762             return false;
     762            goto escape;
    763763        }
    764764
    … …  
    769769                binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);
    770770
    771         float deltaBounds, deltaNums, prevPixels;
    772         float percentNums, base, binLoF32;
    773 
    774         // find the -1 sigma points with linear interpolation
    775         deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
    776         if (binLo == 0) {
    777             deltaNums = cumulative->nums->data.F32[0];
    778             prevPixels = 0;
    779         } else {
    780             deltaNums = cumulative->nums->data.F32[binLo] - cumulative->nums->data.F32[binLo-1];
    781             prevPixels = cumulative->nums->data.F32[binLo-1];
    782         }
    783         percentNums = (totalDataPoints * 0.158655f) - prevPixels;
    784         base = cumulative->bounds->data.F32[binLo];
    785         binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
    786         psTrace("psLib.math", 6,
    787                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
    788                 base, deltaBounds, deltaNums, prevPixels, percentNums);
    789 
    790         // find the +1 sigma points with linear interpolation
    791         deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi];
    792         if (binHi == 0) {
    793             deltaNums = cumulative->nums->data.F32[0];
    794             prevPixels = 0;
    795         } else {
    796             deltaNums = cumulative->nums->data.F32[binHi] - cumulative->nums->data.F32[binHi-1];
    797             prevPixels = cumulative->nums->data.F32[binHi-1];
    798         }
    799         percentNums = (totalDataPoints * 0.841345f) - prevPixels;
    800         base = cumulative->bounds->data.F32[binHi];
    801         float binHiF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 84.1345% mark
    802         psTrace("psLib.math", 6,
    803                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
    804                 base, deltaBounds, deltaNums, prevPixels, percentNums);
     771        // find the +0.5 and -0.5 sigma points with linear interpolation.  binLo and binHi are the bins
     772        // containing the value of interest.  constrain the answer to be within the current bin.
     773        // (extrapolation should not be needed and will result in errors)
     774        float binLoF32, binHiF32, binL2F32, binH2F32;
     775        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f);
     776        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f);
     777        PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f);
     778        PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f);
    805779
    806780        // report +/- 1 sigma points
    … …  
    809783                binLoF32, binHiF32);
    810784
    811         // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions.
    812         sigma = (binHiF32 - binLoF32) / 2.0;
     785        // ADD step 5: Determine SIGMA as (1/2 of) the distance between these positions.
     786        // sigma = (binHiF32 - binLoF32) / 2.0;
     787        sigma = (binH2F32 - binL2F32);
     788
    813789        psTrace("psLib.math", 6, "The current sigma is %f.\n", sigma);
    814790        stats->robustStdev = sigma;
    … …  
    835811            // Free the histograms; they will be recreated on the next iteration, with new bounds
    836812            psFree(histogram);
     813            histogram = NULL;
     814
    837815            psFree(cumulative);
     816            cumulative = NULL;
    838817        } else {
    839818            // We've got the bin size correct now
    … …  
    842821        }
    843822    }
    844     psFree(histogram);
     823
     824    // XXX test lines while studying algorithm errors
     825    // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld  : %f %f %f\n",
     826    // stats->robustMedian, stats->robustStdev,
     827    // binLo, binL2, binMedian, binH2, binHi, binSize, max, min);
    845828
    846829    // ADD step 7: Find the bins which contains the 25% and 75% data points.
    847     long binLo25 = 0;                   // Bin containing the 25% value
    848     long binHi25 = 0;                   // Bin containing the 75% value
    849 
    850     tmpScalar.data.F32 = totalDataPoints * 0.25f;
    851     if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
    852         // Special case where its in first bin.  The last bin should be handled automatically.
    853         binLo25 = 0;
    854     } else {
    855         binLo25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    856     }
    857     tmpScalar.data.F32 = totalDataPoints * 0.75f;
    858     if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
    859         // Special case where its in first bin.  The last bin should be handled automatically.
    860         binHi25 = 0;
    861     } else {
    862         binHi25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    863     }
    864     if ((binLo25 < 0) || (binHi25 < 0)) {
    865         psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n");
    866         psFree(statsMinMax);
    867         psFree(cumulative);
    868         psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    869         psFree(mask);
    870         return false;
    871     }
     830    long binLo25, binHi25;
     831    PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
     832    PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
    872833    psTrace("psLib.math", 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);
    873834
    … …  
    875836    // positions.
    876837    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
    877                       binLo25, totalDataPoints * 0.25f);
     838                                                        binLo25, totalDataPoints * 0.25f);
    878839    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
    879                        binHi25, totalDataPoints * 0.75f);
    880     psFree(cumulative);
    881 
     840                                                         binHi25, totalDataPoints * 0.75f);
    882841    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    883842        psError(PS_ERR_UNKNOWN, false,
    884843                "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
    885         psFree(statsMinMax);
    886         psTrace("psLib.math", 4, "---- %s(1) end  ----\n", __func__);
    887         psFree(mask);
    888         return false;
    889     }
    890 
    891     stats->results |= PS_STAT_ROBUST_QUARTILE;
     844        goto escape;
     845    }
     846
    892847    stats->robustLQ = binLo25F32;
    893848    stats->robustUQ = binHi25F32;
    … …  
    897852    for (long i = 0 ; i < myVector->n ; i++) {
    898853        if (!mask->data.U8[i] &&
    899                 (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
     854            (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
    900855            N50++;
    901856        }
    … …  
    905860
    906861    // Clean up
     862    psFree(histogram);
     863    psFree(cumulative);
    907864    psFree(statsMinMax);
    908865    psFree(mask);
    … …  
    910867    stats->results |= PS_STAT_ROBUST_MEDIAN;
    911868    stats->results |= PS_STAT_ROBUST_STDEV;
     869    stats->results |= PS_STAT_ROBUST_QUARTILE;
    912870
    913871    psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
    914872    return true;
     873
     874escape:
     875    stats->robustMedian = NAN;
     876    stats->robustStdev = NAN;
     877    stats->robustUQ = NAN;
     878    stats->robustLQ = NAN;
     879    stats->robustN50 = 0;
     880    stats->results |= PS_STAT_ROBUST_MEDIAN;
     881    stats->results |= PS_STAT_ROBUST_STDEV;
     882    stats->results |= PS_STAT_ROBUST_QUARTILE;
     883
     884    psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
     885
     886    psFree(histogram);
     887    psFree(cumulative);
     888    psFree(statsMinMax);
     889    psFree(mask);
     890
     891    return false;
    915892}
    916893
    … …  
    943920    float guessMean = stats->robustMedian;  // pass the guess mean
    944921
     922    psVectorBinaryDisectResult result;
     923
    945924    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
    946925    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
    … …  
    950929        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    951930
    952         psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
     931        psScalar tmpScalar;                 // Temporary scalar variable, for psVectorBinaryDisect
    953932        tmpScalar.type.type = PS_TYPE_F32;
    954933
    … …  
    995974        }
    996975
    997         long binMin = 0;                // Low bin to check
    998         long binMax = 0;                // High bin to check
    999 
    1000976        // Fit a Gaussian to the bins in the requested range about the guess mean
    1001977        // min,max overrides clipSigma
    … …  
    1016992        }
    1017993
    1018         tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
    1019         if (tmpScalar.data.F32 < min) {
    1020             binMin = 0;
    1021         } else {
    1022             binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    1023         }
    1024         tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
    1025         if (tmpScalar.data.F32 > max) {
    1026             binMax = histogram->bounds->n - 1;
    1027         } else {
    1028             binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    1029         }
    1030         if (binMin < 0 || binMax < 0) {
    1031             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
    1032             psFree(statsMinMax);
    1033             psFree(histogram);
    1034             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    1035             return false;
    1036         }
     994        // select the min and max bins, saturating on the lower and upper end-points
     995        long binMin, binMax;
     996        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);
     997        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);
    1037998
    1038999        // Generate the variables that will be used in the Gaussian fitting
    … …  
    11421103    float guessMean = stats->robustMedian;  // pass the guess mean
    11431104
     1105    psVectorBinaryDisectResult result;
    11441106    psTrace("psLib.math", 6, "The ** starting ** guess mean  is %f.\n", guessMean);
    11451107    psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
    … …  
    11491111        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    11501112
    1151         psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
     1113        psScalar tmpScalar;                 // Temporary scalar variable, for psVectorBinaryDisect
    11521114        tmpScalar.type.type = PS_TYPE_F32;
    11531115
    … …  
    11941156        }
    11951157
    1196         long binMin = 0;                // Low bin to check
    1197         long binMax = 0;                // High bin to check
    1198 
    11991158        // Fit a Gaussian to the bins in the requested range about the guess mean
    12001159        // min,max overrides clipSigma
    … …  
    12151174        }
    12161175
    1217         tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
    1218         if (tmpScalar.data.F32 < min) {
    1219             binMin = 0;
    1220         } else {
    1221             binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    1222         }
    1223         tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
    1224         if (tmpScalar.data.F32 > max) {
    1225             binMax = histogram->bounds->n - 1;
    1226         } else {
    1227             binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    1228         }
    1229         if (binMin < 0 || binMax < 0) {
    1230             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
    1231             psFree(statsMinMax);
    1232             psFree(histogram);
    1233             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    1234             return false;
    1235         }
     1176        // select the min and max bins, saturating on the lower and upper end-points
     1177        long binMin, binMax;
     1178        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);
     1179        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);
    12361180
    12371181        // Generate the variables that will be used in the Gaussian fitting
    … …  
    13551299    float guessMean = stats->robustMedian;  // pass the guess mean
    13561300
     1301    psVectorBinaryDisectResult result;
    13571302    psTrace("psLib.math", 6, "The ** starting ** guess mean  is %f.\n", guessMean);
    13581303    psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
    … …  
    13621307        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    13631308
    1364         psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
     1309        psScalar tmpScalar;                 // Temporary scalar variable, for psVectorBinaryDisect
    13651310        tmpScalar.type.type = PS_TYPE_F32;
    13661311
    … …  
    14091354        // now fit a Gaussian to the upper and lower halves about the peak independently
    14101355
    1411         long binMin, binMax;
    1412         float upperMean, upperStdev;
    1413 
    14141356        // set the full-range upper and lower limits
    14151357        psF32 maxFitSigma = 2.0;
    … …  
    14291371        }
    14301372
    1431         tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
    1432         if (tmpScalar.data.F32 < min) {
    1433             binMin = 0;
    1434         } else {
    1435             binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    1436         }
    1437         tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
    1438         if (tmpScalar.data.F32 > max) {
    1439             binMax = histogram->bounds->n;
    1440         } else {
    1441             binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar) + 1;
    1442         }
    1443         if (binMin < 0 || binMax < 0) {
    1444             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
    1445             psFree(statsMinMax);
    1446             psFree(histogram);
    1447             psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
    1448             return false;
    1449         }
     1373        // select the min and max bins, saturating on the lower and upper end-points
     1374        long binMin, binMax;
     1375        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);
     1376        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);
    14501377
    14511378        // search for mode (peak of histogram within range mean-2sigma - mean+2sigma
    … …  
    14631390        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", binMax - binMin);
    14641391        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
    1465         psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
     1392        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1);
    14661393        psTrace("psLib.math", 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak);
    14671394
    … …  
    15811508
    15821509            // calculate upper mean & stdev from parabolic fit -- use this as the result
    1583             upperStdev = sqrt(-0.5/poly->coeff[2]);
    1584             upperMean = poly->coeff[1]*PS_SQR(upperStdev);
     1510            float upperStdev = sqrt(-0.5/poly->coeff[2]);
     1511            float upperMean = poly->coeff[1]*PS_SQR(upperStdev);
    15851512            psTrace("psLib.math", 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    15861513            psTrace("psLib.math", 6, "The upper mean  is %f.\n", upperMean);
    … …  
    16801607            // took way too long.
    16811608            //
    1682             #if 0
     1609#if 0
    16831610
    16841611            gaussWidth = 10;
    1685             #endif
     1612#endif
    16861613            //
    16871614            // We determine the bin numbers (jMin:jMax) corresponding to a
    … …  
    19261853psStatsOptions psStatsOptionFromString(const char *string)
    19271854{
    1928     #define READ_STAT(NAME, SYMBOL) \
     1855#define READ_STAT(NAME, SYMBOL) \
    19291856    if (strcasecmp(string, NAME) == 0) { \
    19301857        return SYMBOL; \
    … …  
    19641891    psString string = NULL;             // String to return
    19651892
    1966     #define WRITE_STAT(NAME, SYMBOL) \
     1893#define WRITE_STAT(NAME, SYMBOL) \
    19671894    if (option & SYMBOL) { \
    19681895        psStringAppend(&string, "%s ", NAME); \
    … …  
    20231950{
    20241951    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
    2025     case PS_STAT_SAMPLE_MEAN:
    2026     case PS_STAT_SAMPLE_MEDIAN:
    2027     case PS_STAT_SAMPLE_STDEV:
    2028     case PS_STAT_SAMPLE_QUARTILE:
    2029     case PS_STAT_ROBUST_MEDIAN:
    2030     case PS_STAT_ROBUST_STDEV:
    2031     case PS_STAT_ROBUST_QUARTILE:
    2032     case PS_STAT_FITTED_MEAN:
    2033     case PS_STAT_FITTED_STDEV:
    2034     case PS_STAT_FITTED_MEAN_V2:
    2035     case PS_STAT_FITTED_STDEV_V2:
    2036     case PS_STAT_FITTED_MEAN_V3:
    2037     case PS_STAT_FITTED_STDEV_V3:
    2038     case PS_STAT_CLIPPED_MEAN:
    2039     case PS_STAT_CLIPPED_STDEV:
    2040     case PS_STAT_MAX:
    2041     case PS_STAT_MIN:
     1952      case PS_STAT_SAMPLE_MEAN:
     1953      case PS_STAT_SAMPLE_MEDIAN:
     1954      case PS_STAT_SAMPLE_STDEV:
     1955      case PS_STAT_SAMPLE_QUARTILE:
     1956      case PS_STAT_ROBUST_MEDIAN:
     1957      case PS_STAT_ROBUST_STDEV:
     1958      case PS_STAT_ROBUST_QUARTILE:
     1959      case PS_STAT_FITTED_MEAN:
     1960      case PS_STAT_FITTED_STDEV:
     1961      case PS_STAT_FITTED_MEAN_V2:
     1962      case PS_STAT_FITTED_STDEV_V2:
     1963      case PS_STAT_FITTED_MEAN_V3:
     1964      case PS_STAT_FITTED_STDEV_V3:
     1965      case PS_STAT_CLIPPED_MEAN:
     1966      case PS_STAT_CLIPPED_STDEV:
     1967      case PS_STAT_MAX:
     1968      case PS_STAT_MIN:
    20421969        return option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE);
    2043     default:
     1970      default:
    20441971        return 0;
    20451972    }
    … …  
    20521979    // We could call psStatsSingle to check, but it would be a waste since we effectively do it anyway
    20531980    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
    2054     case PS_STAT_SAMPLE_MEAN:
     1981      case PS_STAT_SAMPLE_MEAN:
    20551982        return stats->sampleMean;
    2056     case PS_STAT_SAMPLE_MEDIAN:
     1983      case PS_STAT_SAMPLE_MEDIAN:
    20571984        return stats->sampleMedian;
    2058     case PS_STAT_SAMPLE_STDEV:
     1985      case PS_STAT_SAMPLE_STDEV:
    20591986        return stats->sampleStdev;
    2060     case PS_STAT_ROBUST_MEDIAN:
     1987      case PS_STAT_ROBUST_MEDIAN:
    20611988        return stats->robustMedian;
    2062     case PS_STAT_ROBUST_STDEV:
     1989      case PS_STAT_ROBUST_STDEV:
    20631990        return stats->robustStdev;
    2064     case PS_STAT_FITTED_MEAN:
     1991      case PS_STAT_FITTED_MEAN:
    20651992        return stats->fittedMean;
    2066     case PS_STAT_FITTED_STDEV:
     1993      case PS_STAT_FITTED_STDEV:
    20671994        return stats->fittedStdev;
    2068     case PS_STAT_FITTED_MEAN_V2:
     1995      case PS_STAT_FITTED_MEAN_V2:
    20691996        return stats->fittedMean;
    2070     case PS_STAT_FITTED_STDEV_V2:
     1997      case PS_STAT_FITTED_STDEV_V2:
    20711998        return stats->fittedStdev;
    2072     case PS_STAT_FITTED_MEAN_V3:
     1999      case PS_STAT_FITTED_MEAN_V3:
    20732000        return stats->fittedMean;
    2074     case PS_STAT_FITTED_STDEV_V3:
     2001      case PS_STAT_FITTED_STDEV_V3:
    20752002        return stats->fittedStdev;
    2076     case PS_STAT_CLIPPED_MEAN:
     2003      case PS_STAT_CLIPPED_MEAN:
    20772004        return stats->clippedMean;
    2078     case PS_STAT_CLIPPED_STDEV:
     2005      case PS_STAT_CLIPPED_STDEV:
    20792006        return stats->clippedStdev;
    2080     case PS_STAT_MAX:
     2007      case PS_STAT_MAX:
    20812008        return stats->max;
    2082     case PS_STAT_MIN:
     2009      case PS_STAT_MIN:
    20832010        return stats->min;
    2084     case PS_STAT_SAMPLE_QUARTILE:
    2085     case PS_STAT_ROBUST_QUARTILE:
     2011      case PS_STAT_SAMPLE_QUARTILE:
     2012      case PS_STAT_ROBUST_QUARTILE:
    20862013        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot return a single quartile value; "
    20872014                "get them yourself.\n");
    20882015        return NAN;
    2089     default:
     2016      default:
    20902017        return NAN;
    20912018    }
    … …  
    21022029                              psF32 xLo,
    21032030                              psF32 xHi
    2104                              )
     2031    )
    21052032{
    21062033    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
    … …  
    21302057*****************************************************************************/
    21312058static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
    2132         psVector *yVec,
    2133         psS32 binNum,
    2134         psF32 yVal
    2135                                               )
     2059                                               psVector *yVec,
     2060                                               psS32 binNum,
     2061                                               psF32 yVal
     2062    )
    21362063{
    21372064    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    21712098        //
    21722099        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
    2173                 ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2100               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
    21742101            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    21752102                    _("Specified yVal, %g, is not within y-range, %g to %g."),
    … …  
    21812108        //
    21822109        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
    2183                 ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
     2110            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
    21842111            psError(PS_ERR_UNKNOWN, true,
    21852112                    "This routine must be called with monotonically increasing or decreasing data points.\n");
    … …  
    22522179                                   const psVector *params,
    22532180                                   const psVector *coords
    2254                                   )
     2181    )
    22552182{
    22562183    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
Note: See TracChangeset for help on using the changeset viewer.