IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 1, 2013, 1:59:06 PM (13 years ago)
Author:
bills
Message:

Revert psStats.c to r35815, the version in the ipp-20130712 production tag
to fix problems with background modeling in psphot.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c

    r36256 r36262  
    140140}
    141141
    142 // Debug information
    143 #define CZW 0
    144142
    145143/*****************************************************************************/
    … …  
    174172*****************************************************************************/
    175173
    176 // static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    177 static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     174static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    178175
    179176/******************************************************************************
    … …  
    232229        }
    233230        count++;
    234 
    235231    }
    236232    if (errors) {
    … …  
    797793        } else {
    798794            // Determine the bin size of the robust histogram, using the pre-defined number of bins
    799             binSize = (max - min) / INITIAL_NUM_BINS;
     795            binSize = (max - min) / INITIAL_NUM_BINS;
    800796        }
    801797        psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
    … …  
    880876        cumulative = psHistogramAlloc(min, max, numBins);
    881877        cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
    882         cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1];
    883 
    884         // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
    885         long Nc = 1;  // track the current bin of cumulative
    886         // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the
    887         // upper end of the current histogram bin
    888         for (long i = 1; i < histogram->nums->n - 1; i++) {
    889             if (histogram->nums->data.F32[i] == 0.0) continue;
    890             cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
    891             cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1];
    892             Nc ++;
    893         }
    894         long Nlast = Nc - 1;  // last valid cumulative bin
    895         for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
    896             cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
    897             cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
    898         }
    899        
     878        for (long i = 1; i < histogram->nums->n; i++) {
     879            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
     880            cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
     881        }
    900882        if (psTraceGetLevel("psLib.math") >= 8) {
    901883            PS_VECTOR_PRINT_F32(cumulative->bounds);
    … …  
    913895
    914896        // ADD step 3: Interpolate to the exact 50% position in bin units
    915         // stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    916         // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     897        stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     898        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    917899        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
    918         // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there.
    919         // Instead, do a 5-point linear fit.
    920         stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    921900
    922901        // convert bin to bin value: this is the robust histogram median.
    … …  
    933912        PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0);
    934913        PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0);
    935         PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0);
    936         PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0);
    937        
    938        
     914        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022481f, 0);
     915        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0);
     916
    939917        psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
    940918                binLo, binHi);
    … …  
    948926            goto escape;
    949927        }
    950    
     928
    951929        // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma
    952930        // positions.
    … …  
    960938        // (extrapolation should not be needed and will result in errors)
    961939        float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
    962 #if (0)
    963940        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo,
    964941                            totalDataPoints * 0.158655f);
    … …  
    970947                            totalDataPoints * 0.691462f);
    971948        PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4,
    972                             totalDataPoints * 0.022750f);
     949                            totalDataPoints * 0.022481f);
    973950        PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4,
    974                             totalDataPoints * 0.977250f);
    975 #else
    976         binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655);
    977         binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345);         
    978         binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538);
    979         binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462);         
    980         binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750);
    981         binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250);
    982 #endif 
     951                            totalDataPoints * 0.977519f);
     952
    983953        // report +/- 1 sigma points
    984954        psTrace(TRACE, 5,
    … …  
    989959                binL2F32, binH2F32);
    990960        psTrace(TRACE, 5,
    991                 "The exact 02.2275 and 97.7250 percent data point positions are: (%f, %f)\n",
     961                "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n",
    992962                binL4F32, binH4F32);
    993963
    994         // If some of the fits failed, attempt to fix this
    995         if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; }
    996         if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; }
    997         if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; }
    998         if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; }
    999         if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; }
    1000         if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; }
    1001        
    1002964        // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma)
    1003 
    1004 
    1005965        float sigma1 = (binH2F32 - binL2F32);
    1006966        float sigma2 = (binHiF32 - binLoF32) / 2.0;
    1007967        float sigma4 = (binH4F32 - binL4F32) / 4.0;
    1008968
    1009         // Fix again?
    1010         if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; }
    1011         if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; }
    1012         if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.0; }
    1013        
    1014969        // take the smallest of the three: if we have a clump with wide outliers, sigma2 and
    1015970        // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2
    1016971        // will be biased high.
    1017         //        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
    1018         // CZW: Instead, take the median.  Taking the MIN forces a bias on unbiased data.
    1019         //      It seems like occasionally getting the wrong answer on a complex distribution
    1020         //      is more acceptable than always getting the wrong answer for simple ones.
    1021 
    1022        
    1023         sigma = PS_MAX( PS_MIN(sigma1,sigma2),
    1024                         PS_MIN( PS_MAX(sigma1,sigma2),
    1025                                 sigma4));
     972        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
    1026973
    1027974        psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1);
    … …  
    1030977
    1031978        psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
    1032         //        stats->robustStdev = sigma;
    1033         stats->robustStdev = sigma;
    1034 
    1035 #if (CZW && 0)
    1036         // Skewness check: Find least biased sample for each pair.
    1037         sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian,
    1038                               stats->robustMedian - binL2F32);
    1039         sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian,
    1040                               stats->robustMedian - binLoF32);
    1041         sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian,
    1042                               stats->robustMedian - binL4F32);
    1043         // Kurtosis check: Take median sample as the solution.
    1044         stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2),
    1045                                      PS_MIN( PS_MAX(sigma1,sigma2),
    1046                                              sigma4));
    1047 #endif
    1048 
    1049        
    1050 #if (CZW)
    1051         printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
    1052                binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
    1053                sigma1,sigma2,sigma4,sigma);
    1054        
    1055         printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
    1056                iterate,
    1057            stats->robustMedian,stats->robustStdev,
    1058            fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
    1059 
    1060            cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
    1061            cumulative->bounds->data.F32[binMedian-1],
    1062            cumulative->bounds->data.F32[binMedian],
    1063            cumulative->bounds->data.F32[binMedian+1],
    1064            cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
    1065 
    1066            cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
    1067            cumulative->nums->data.F32[binMedian-1],
    1068            cumulative->nums->data.F32[binMedian],
    1069            cumulative->nums->data.F32[binMedian+1],
    1070            cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
    1071         //      PS_VECTOR_PRINT_F32(histogram->bounds);
    1072         //      PS_VECTOR_PRINT_F32(histogram->nums);
    1073         //      PS_VECTOR_PRINT_F32(cumulative->bounds);
    1074         //      PS_VECTOR_PRINT_F32(cumulative->nums);
    1075 #endif
     979        stats->robustStdev = sigma;
    1076980
    1077981        // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
    … …  
    1079983        if (sigma < (3.0 * binSize)) {
    1080984            psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
    1081 
    1082             // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
    1083             psF32 medianLo = stats->robustMedian - 25*binSize;
    1084             psF32 medianHi = stats->robustMedian + 25*binSize;
    1085 
    1086             // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
    1087             // long maskHi = PS_MIN(cumulative->bounds->n - 1, (binMedian + 25)); // High index for masking
    1088             // psF32 medianLo = cumulative->bounds->data.F32[maskLo]; // Value at low index
    1089             // psF32 medianHi = cumulative->bounds->data.F32[maskHi]; // Value at high index
     985            long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
     986            long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
     987            psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
     988            psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
    1090989            psTrace(TRACE, 6, "Masking data more than 25 bins from the median\n");
    1091             // psTrace(TRACE, 6, "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi);
     990            psTrace(TRACE, 6,
     991                    "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n",
     992                    binMedian, maskLo, maskHi);
    1092993            psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
    1093             int Nmasked = 0;
    1094994            for (long i = 0 ; i < myVector->n ; i++) {
    1095995                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
    1096                     if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue;
     996                    // XXXX is this correct?  is MASK_MARK safe?
    1097997                    mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
    1098998                    psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
    1099                     Nmasked ++;
    1100999                }
    11011000            }
    1102 
    1103             if (Nmasked == 0) {
    1104                 // no significant change to the sigma & binsize -- we are done here
    1105                 iterate = -1;
    1106                 continue;
    1107             }
    11081001
    11091002            // Free the histograms; they will be recreated on the next iteration, with new bounds
    … …  
    11371030        }
    11381031    }
    1139    
     1032
    11401033    // XXX test lines while studying algorithm errors
    11411034    // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld  : %f %f %f\n",
    … …  
    11471040    PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
    11481041    PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
    1149     psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25);
     1042    psTrace(TRACE, 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);
    11501043
    11511044    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
    11521045    // positions.
    1153     psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
    1154     psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
     1046    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
     1047    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
    11551048    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    1156         COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n");
     1049        COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
    11571050        goto escape;
    11581051    }
    … …  
    12071100 * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4
    12081101 * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev
    1209  * robustN50 should also be set gaussian fit is performed using 1D polynomial to ln(y) this
     1102 * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this
    12101103 * version follows the upper portion of the distribution until it passes 0.5*peak
    12111104 ********************/
    … …  
    12421135        return true;
    12431136    }
    1244     if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); }
     1137
    12451138    float guessStdev = stats->robustStdev;  // pass the guess sigma
    12461139    float guessMean = stats->robustMedian;  // pass the guess mean
    … …  
    12621155            // set roughly so that the lowest bins have about 2 cnts
    12631156            // Nsmallest ~ N50 / (4*dN))
    1264             psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); 
    1265            binSize = guessStdev / dN;
     1157            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
     1158            binSize = guessStdev / dN;
    12661159        }
    12671160
    … …  
    12891182        // XXX can we calculate the binMin, binMax **before** building this histogram?
    12901183        // if the range is too absurd, adjust numBins & binSize
    1291         // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
    1292         // a number of bins, and then set the min/max values to put those bins sanely around the mean.
    12931184        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
    1294         //        binSize = (max - min) / (float) numBins;
     1185        binSize = (max - min) / (float) numBins;
    12951186        psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
    12961187        psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
    12971188        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    12981189
    1299 
    1300 #define FITTED_CLIPPING_NUM 5.0
    1301         if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) {
    1302           min = guessMean - FITTED_CLIPPING_NUM * guessStdev;
    1303         }
    1304         if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) {
    1305           max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
    1306         }
    1307        
    13081190        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    13091191        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    … …  
    13401222        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
    13411223        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    1342 
    13431224        if (binMin == binMax) {
    13441225            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    … …  
    13671248        psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
    13681249
    1369        
    13701250        float lowfitMean = NAN;
    13711251        float lowfitStdev = NAN;
    … …  
    14051285            }
    14061286            y->n = x->n = j;
    1407            
     1287
    14081288            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    14091289            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
    … …  
    14171297                psErrorClear();
    14181298                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
    1419 
    14201299                psFree(poly);
    14211300                psFree(histogram);
    … …  
    14991378            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    15001379            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    1501 #if (CZW && 0)
    1502             for (long i = 0; i < x->n; i++) {
    1503               printf("CZW: Dcheck: %ld %f %f %f\n",
    1504                      i,x->data.F32[i],y->data.F32[i],
    1505                      poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
    1506                      poly->coeff[2] * pow(x->data.F32[i],2));
    1507             }
    1508 #endif
    15091380            psFree(x);
    15101381            psFree(y);
    … …  
    15221393            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
    15231394            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
    1524 
    15251395#ifndef PS_NO_TRACE
    15261396            psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    … …  
    15451415            }
    15461416
    1547            
    15481417            psFree (poly);
    15491418        }
    … …  
    15681437            done = true;
    15691438        }
    1570 
    1571        
    1572 #if (CZW && 1)
    1573         printf("CZW IN FITTED: iter   %d %f \n"
    1574                "               low    %f %f \n"
    1575                "               full   %f %f \n"
    1576                "               robust %f %f \n"
    1577                "               final  %f %f\n",
    1578                iteration,minValueSym,
    1579                lowfitMean,lowfitStdev,
    1580                fullfitMean,fullfitStdev,
    1581                stats->robustMedian,stats->robustStdev,
    1582                guessMean,guessStdev);
    1583 #endif
    15841439
    15851440        // Clean up after fitting
    … …  
    21101965// other private functions used above
    21111966
    2112 # if (0)
    21131967static psF32 QuadraticInverse(psF32 a,
    21141968                              psF32 b,
    … …  
    21321986    return 0.5 * (xLo + xHi);
    21331987}
    2134 
    2135 static psF32 LinearInverse(psF32 a,
    2136                            psF32 b,
    2137                            psF32 y,
    2138                            psF32 xLo,
    2139                            psF32 xHi
    2140     )
    2141 {
    2142     psF64 x = (y - b) / a;
    2143 
    2144     if (xLo <= x && x <= xHi) {
    2145         return x;
    2146     }
    2147     return 0.5 * (xLo + xHi);
    2148 }
    2149 # endif
    21501988
    21511989# if (0)
    … …  
    24042242    return tmpFloat;
    24052243}
     2244# endif
    24062245
    24072246/******************************************************************************
    … …  
    24372276    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
    24382277
    2439     //    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    2440     //    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    2441     psVector *x = psVectorAlloc(5, PS_TYPE_F64);
    2442     psVector *y = psVectorAlloc(5, PS_TYPE_F64);
     2278    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     2279    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    24432280    psF32 tmpFloat = 0.0f;
    24442281
    2445     //    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
    2446     if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
     2282    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
    24472283        // The general case.  We have all three points.
    2448       //        x->data.F64[0] = binNum - 1;
    2449       //        x->data.F64[1] = binNum;
    2450       //        x->data.F64[2] = binNum + 1;
    2451       x->data.F64[0] = xVec->data.F32[binNum - 2];
    2452       x->data.F64[1] = xVec->data.F32[binNum - 1];
    2453       x->data.F64[2] = xVec->data.F32[binNum + 0];
    2454       x->data.F64[3] = xVec->data.F32[binNum + 1];
    2455       x->data.F64[4] = xVec->data.F32[binNum + 2];
    2456         y->data.F64[0] = yVec->data.F32[binNum - 2];
    2457         y->data.F64[1] = yVec->data.F32[binNum - 1];
    2458         y->data.F64[2] = yVec->data.F32[binNum + 0];
    2459         y->data.F64[3] = yVec->data.F32[binNum + 1];
    2460         y->data.F64[4] = yVec->data.F32[binNum + 2];
     2284        x->data.F64[0] = binNum - 1;
     2285        x->data.F64[1] = binNum;
     2286        x->data.F64[2] = binNum + 1;
     2287        y->data.F64[0] = yVec->data.F32[binNum - 1];
     2288        y->data.F64[1] = yVec->data.F32[binNum];
     2289        y->data.F64[2] = yVec->data.F32[binNum + 1];
    24612290        psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
    24622291        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
    24632292        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    24642293
    2465 
    24662294        // Ensure that the y value lies within range of the y values.
    2467         if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
    2468                ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2295        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
     2296               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
    24692297            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    24702298                    _("Specified yVal, %g, is not within y-range, %g to %g."),
    … …  
    25032331
    25042332        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
    2505         float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]);
     2333        float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
    25062334        psFree(myPoly);
    25072335
    … …  
    25132341            return(NAN);
    25142342        }
    2515        
     2343
    25162344        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
    2517         //      assert (binValue >= binNum - 1);
    2518         //      assert (binValue <= binNum + 1);
    2519 
    2520         //      int fitBin = binValue;
    2521         //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
    2522         //        float dY = binValue - fitBin;
    2523         //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
    2524         tmpFloat = binValue;
    2525        
     2345        assert (binValue >= binNum - 1);
     2346        assert (binValue <= binNum + 1);
     2347
     2348        int fitBin = binValue;
     2349        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2350        float dY = binValue - fitBin;
     2351        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
    25262352    } else {
    25272353        // These are special cases where the bin is at the beginning or end of the vector.
    … …  
    25552381    return tmpFloat;
    25562382}
    2557 # endif
    2558 
    2559 
    2560 /******************************************************************************
    2561 fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
    2562 which fits a quadratic to three points and returns the x bin value corresponding to the input
    2563 y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
    2564 points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
    2565 x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
    2566 determines for what value x does that quadratic f(x) = yVal (the input parameter).
    2567 
    2568 XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation
    2569 could easily be done with statically allocated doubles, skipping the psLib versions of
    2570 polynomial fitting, etc.
    2571 
    2572 *****************************************************************************/
    2573 static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec,
    2574                                               psVector *yVec,
    2575                                               psS32 binNum,
    2576                                               psF32 yVal
    2577     )
    2578 {
    2579 
    2580 # if (1)
    2581 # define HALF_SIZE 2
    2582   double Sx = 0.0;
    2583 
    2584   double Sy = 0.0;
    2585   double Sxx = 0.0;
    2586   double Sxy = 0.0;
    2587   double deltaY = 0.0;
    2588   int N = 0;
    2589 
    2590   for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
    2591     if ((u >= 0)&&(u < yVec->n)) {
    2592       if (u+1 < xVec->n) {
    2593         Sx += yVec->data.F32[u];
    2594         Sxx += PS_SQR(yVec->data.F32[u]);
    2595 
    2596         deltaY = xVec->data.F32[u];
    2597         //deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
    2598         Sy += deltaY;
    2599         Sxy += yVec->data.F32[u] * deltaY;
    2600         N += 1;
    2601       }
    2602     }
    2603   }
    2604   double Det = N * Sxx - Sx * Sx;
    2605   if (Det == 0.0) return NAN;
    2606   if (N == 0) return NAN;
    2607 
    2608   double C0 = (Sy*Sxx - Sx*Sxy) / Det;
    2609   double C1 = (Sxy*N - Sx*Sy) / Det;
    2610  
    2611   double value = C0 + yVal*C1;
    2612   return value;
    2613  
    2614  
    2615 # else
    2616     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
    2617     if (psTraceGetLevel("psLib.math") >= 8) {
    2618         PS_VECTOR_PRINT_F32(xVec);
    2619         PS_VECTOR_PRINT_F32(yVec);
    2620     }
    2621 
    2622     PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
    2623     PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
    2624     PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
    2625     PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
    2626     PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
    2627     PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
    2628 
    2629     //    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    2630     //    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    2631     psVector *x = psVectorAlloc(5, PS_TYPE_F64);
    2632     psVector *y = psVectorAlloc(5, PS_TYPE_F64);
    2633     psF32 tmpFloat = 0.0f;
    2634 
    2635     if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
    2636         x->data.F64[0] = xVec->data.F32[binNum - 2];
    2637         x->data.F64[1] = xVec->data.F32[binNum - 1];
    2638         x->data.F64[2] = xVec->data.F32[binNum + 0];
    2639         x->data.F64[3] = xVec->data.F32[binNum + 1];
    2640         x->data.F64[4] = xVec->data.F32[binNum + 2];
    2641 
    2642         y->data.F64[0] = yVec->data.F32[binNum - 2];
    2643         y->data.F64[1] = yVec->data.F32[binNum - 1];
    2644         y->data.F64[2] = yVec->data.F32[binNum + 0];
    2645         y->data.F64[3] = yVec->data.F32[binNum + 1];
    2646         y->data.F64[4] = yVec->data.F32[binNum + 2];
    2647         psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
    2648         psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
    2649         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    2650 
    2651         // Ensure that the y value lies within range of the y values.
    2652         if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
    2653                ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
    2654             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    2655                     _("Specified yVal, %g, is not within y-range, %g to %g."),
    2656                     (psF64)yVal, y->data.F64[0], y->data.F64[2]);
    2657             return NAN;
    2658         }
    2659 
    2660         // Ensure that the y values are monotonic.
    2661         if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
    2662             ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
    2663             psError(PS_ERR_UNKNOWN, true,
    2664                     "This routine must be called with monotonically increasing or decreasing data points.\n");
    2665             psFree(x);
    2666             psFree(y);
    2667             return NAN;
    2668         }
    2669 
    2670         // Determine the coefficients of the polynomial.
    2671         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
    2672         if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
    2673             psError(PS_ERR_UNEXPECTED_NULL, false,
    2674                     _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
    2675                       "Returning NAN."));
    2676             psFree(x);
    2677             psFree(y);
    2678             return NAN;
    2679         }
    2680 
    2681         psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
    2682         psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
    2683         psTrace(TRACE, 6, "Fitted y vec is (%f %f)\n",
    2684                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
    2685                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]));
    2686 
    2687         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
    2688         float binValue = LinearInverse(myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]);
    2689         psFree(myPoly);
    2690 
    2691         if (isnan(binValue)) {
    2692             psError(PS_ERR_UNEXPECTED_NULL,
    2693                     false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
    2694             psFree(x);
    2695             psFree(y);
    2696             return(NAN);
    2697         }
    2698        
    2699         // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
    2700         //      assert (binValue >= binNum - 1);
    2701         //      assert (binValue <= binNum + 1);
    2702 
    2703         //      int fitBin = binValue;
    2704         //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
    2705         //        float dY = binValue - fitBin;
    2706         //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
    2707         tmpFloat = binValue;
    2708                
    2709        
    2710     } else {
    2711         // These are special cases where the bin is at the beginning or end of the vector.
    2712         if (binNum == 0) {
    2713             // We have two points only at the beginning of the vectors x and y.
    2714             // X = (dX/dY)(Y - Yo) + Xo
    2715             float dX = xVec->data.F32[1] - xVec->data.F32[0];
    2716             float dY = yVec->data.F32[1] - yVec->data.F32[0];
    2717             if (dY == 0.0) {
    2718                 tmpFloat = xVec->data.F32[0];
    2719             } else {
    2720                 tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
    2721             }
    2722         } else if (binNum == (xVec->n - 1)) {
    2723             // We have two points only at the end of the vectors x and y.
    2724             // X = (dX/dY)(Y - Yo) + Xo
    2725             float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
    2726             float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
    2727             if (dY == 0.0) {
    2728                 tmpFloat = xVec->data.F32[binNum-1];
    2729             } else {
    2730                 tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
    2731             }
    2732         }
    2733     }
    2734 
    2735     psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
    2736     psFree(x);
    2737     psFree(y);
    2738 
    2739     return tmpFloat;
    2740 # endif
    2741 }
Note: See TracChangeset for help on using the changeset viewer.