IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6315


Ignore:
Timestamp:
Feb 2, 2006, 3:30:14 PM (20 years ago)
Author:
gusciora
Message:

Cleaned and modified the robust stats.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6305 r6315  
    1616 * use ->min and ->max (PS_STAT_USE_RANGE)
    1717 *
    18  *  @version $Revision: 1.165 $ $Name: not supported by cvs2svn $
    19  *  @date $Date: 2006-02-02 21:09:08 $
     18 *  @version $Revision: 1.166 $ $Name: not supported by cvs2svn $
     19 *  @date $Date: 2006-02-03 01:30:14 $
    2020 *
    2121 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    878878    psF32 sumSquares = 0.0;     // temporary variable
    879879    psF32 sumDiffs = 0.0;       // temporary variable
    880     //    psF32 sum1;
    881     //    psF32 sum2;
    882880    psF32 errorDivisor = 0.0f;
    883881
    … …  
    913911        } else {
    914912            for (i = 0; i < myVector->n; i++) {
    915                 if ((stats->min <= myVector->data.F32[i]) &&
    916                         (myVector->data.F32[i] <= stats->max)) {
     913                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
    917914                    diff = myVector->data.F32[i] - mean;
    918915                    sumSquares += (diff * diff);
    … …  
    965962        // data ranges are used correctly.
    966963        if (errors != NULL) {
    967             stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
     964            if (0) {
     965                stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
     966            } else {
     967                countFloat = (psF32)countInt;
     968                stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
     969            }
    968970        } else {
    969971            countFloat = (psF32)countInt;
    … …  
    9981000{
    9991001    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     1002    psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__));
    10001003    psF32 clippedMean = 0.0;    // self-explanatory
    10011004    psF32 clippedStdev = 0.0;   // self-explanatory
    … …  
    10711074        return(-2);
    10721075    }
     1076    psTrace(__func__, 6, "The initial sample median is %f\n", statsTmp->sampleMedian);
    10731077
    10741078    // 2. Compute the sample standard deviation.
    … …  
    10811085        return(-2);
    10821086    }
     1087    psTrace(__func__, 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev);
    10831088
    10841089    // 3. Use the sample median as the first estimator of the mean X.
    … …  
    10901095    // 5. Repeat N (stats->clipIter) times:
    10911096    for (psS32 iter = 0; iter < stats->clipIter; iter++) {
     1097        psTrace(__func__, 6, "------------ Iteration %d ------------\n", iter);
    10921098        // a) Exclude all values x_i for which |x_i - x| > K * stdev
    10931099        if (errors != NULL) {
    10941100            for (psS32 j = 0; j < myVector->n; j++) {
    1095                 if (fabs(myVector->data.F32[j] - clippedMean) >
    1096                         (stats->clipSigma * errors->data.F32[j])) {
     1101                if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) {
    10971102                    tmpMask->data.U8[j] = 0xff;
    10981103                }
    … …  
    11101115        p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp);
    11111116        p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp);
     1117        psTrace(__func__, 6, "The new sample mean is %f\n", statsTmp->sampleMean);
     1118        psTrace(__func__, 6, "The new sample stdev is %f\n", statsTmp->sampleStdev);
    11121119
    11131120        // If the new mean and stdev are NAN, we must exit the loop.
    … …  
    11281135    if (stats->options & PS_STAT_CLIPPED_MEAN) {
    11291136        stats->clippedMean = clippedMean;
     1137        psTrace(__func__, 6, "The final clipped mean is %f\n", clippedMean);
    11301138    }
    11311139    // 8. The last calcuated value of stdev is the cliped stdev.
    11321140    if (stats->options & PS_STAT_CLIPPED_STDEV) {
    11331141        stats->clippedStdev = clippedStdev;
     1142        psTrace(__func__, 6, "The final clipped stdev is %f\n", clippedStdev);
    11341143    }
    11351144
    … …  
    11391148}
    11401149
    1141 
    1142 
    1143 
    1144 /******************************************************************************
    1145 p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine
    1146 takes as input a 1-D polynomial of arbitrary order and a range of x-values for
    1147 which it is defined:  [rangeLow, rangeHigh].  It determines the x-value of
    1148 that polynomial such that f(x) == getThisValue.  This function uses a
    1149 binary-search algorithm on the range and assumes that the polynomial is
    1150 monotonically increasing or decreasing within that range.
    1151  
    1152 XXX: Terminate when f(x)-getThisValue is within some error tolerance.
    1153  
    1154 XXX: Create a 2nd-order polynomial version and solve for X analytically.
    1155  *****************************************************************************/
    1156 psF32 p_ps1DPolyMedian(
    1157     psPolynomial1D* myPoly,
    1158     psF32 rangeLow,
    1159     psF32 rangeHigh,
    1160     psF32 getThisValue)
    1161 {
    1162     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    1163     psTrace(__func__, 4, "range (%f, %f) get (%f)\n", rangeLow, rangeHigh, getThisValue);
    1164     PS_ASSERT_POLY_NON_NULL(myPoly, NAN);
    1165     PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN);
    1166     // We ensure that the requested f(y) value, which is getThisValue, is
    1167     // falls within the range of y-values of the polynomial "myPoly" in the
    1168     // specified x-range (rangeLow:rangeHigh).
    1169     psF32 fLo = psPolynomial1DEval(myPoly, rangeLow);
    1170     psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh);
    1171     psTrace(__func__, 4, "function at endpoints are (%f, %f) get (%f)\n", fLo, fHi, getThisValue);
    1172     if (!((fLo <= getThisValue) && (fHi >= getThisValue))) {
    1173         psError(PS_ERR_UNKNOWN, true, "The requested y value (%f) does not fall within the specified range (%f to %f)\n", getThisValue, fLo, fHi);
    1174         psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
    1175         return(NAN);
    1176     }
    1177 
    1178     psS32 numIterations = 0;
    1179     psF32 midpoint = 0.0;
    1180     psF32 oldMidpoint = 1.0;
    1181     psF32 f = 0.0;
    1182 
    1183     while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
    1184         midpoint = (rangeHigh + rangeLow) / 2.0;
    1185         if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
    1186             return (midpoint);
    1187         }
    1188         oldMidpoint = midpoint;
    1189 
    1190         f = psPolynomial1DEval(myPoly, midpoint);
    1191         printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f);
    1192         if (fabs(f - getThisValue) <= FLT_EPSILON) {
    1193             return (midpoint);
    1194         }
    1195 
    1196         if (f > getThisValue) {
    1197             rangeHigh = midpoint;
    1198         } else {
    1199             rangeLow = midpoint;
    1200         }
    1201         numIterations++;
    1202     }
    1203     psTrace(__func__, 4, "---- %s(%f) end ----\n", __func__, midpoint);
    1204     return (midpoint);
    1205 }
    1206 
    1207 
    1208 
    1209 
    1210 
    1211 
    1212 
    1213 
    1214 
    1215 
    1216 
    1217 
    1218 
    1219 
    1220 
    1221 
    1222 
    1223 
    1224 
    1225 
    1226 
    1227 
    1228 
    1229 
    1230 
    1231 
    1232 
    1233 
    1234 
    1235 
    1236 
    1237 
    1238 
    1239 
    1240 
    1241 
    1242 
    1243 
    1244 #define PS_SQRT(x) sqrt(x)
    12451150static psF32 QuadraticInverse(
    12461151    psF32 a,
    12471152    psF32 b,
    12481153    psF32 c,
    1249     psF32 y)
     1154    psF32 y,
     1155    psF32 xLo,
     1156    psF32 xHi)
    12501157{
    12511158    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
    … …  
    12541161    psF64 x2 = -b/(2.0*a) - tmp;
    12551162
    1256     psF64 y1 = (a * x1 * x1) + (b * x1) + c;
    1257     psF64 y2 = (a * x2 * x2) + (b * x2) + c;
    1258 
    1259     printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
    1260     printf("QuadraticInverse: y is %f\n", y);
    1261     printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
    1262     printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
    1263 
    1264     return(x1);
     1163    if (0) {
     1164        psF64 y1 = (a * x1 * x1) + (b * x1) + c;
     1165        psF64 y2 = (a * x2 * x2) + (b * x2) + c;
     1166        printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
     1167        printf("QuadraticInverse: y is %f\n", y);
     1168        printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
     1169        printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
     1170    }
     1171
     1172    if ((xLo <= x1) && (x1 <= xHi)) {
     1173        return(x1);
     1174    } else if ((xLo <= x2) && (x2 <= xHi)) {
     1175        return(x2);
     1176    } else {
     1177        return(0.5 * (xLo + xHi));
     1178    }
    12651179}
    12661180
    … …  
    13891303                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
    13901304
    1391         // Call p_ps1DPolyMedian(), which does a binary search on the
    1392         // polynomial, looking for the value x such that f(x) = yVal
    13931305        psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
    1394         printf("Calling p_ps1DPolyMedian() for the y-value (%.2f) that has an x in the range (%.2f - %.2f)\n", yVal, x->data.F64[0], x->data.F64[2]);
    1395         tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
    1396         printf("Cool, tmpFloat is %.2f\n", tmpFloat);
    1397 
    1398         QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal);
     1306        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
    13991307        psFree(myPoly);
    1400 
    14011308
    14021309        if (isnan(tmpFloat)) {
    … …  
    14701377    deriv->data.F32[1] = tmp / (stdev * stdev * stdev);
    14711378
    1472     // printf("f(x, mean, stdev) = f(%.2f %.2f %.2f) is %.2f\n", x, mean, stdev, psGaussian(x, mean, stdev, false));
    1473 
    14741379    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    14751380    return(psGaussian(x, mean, stdev, false));
    14761381}
    1477 
    1478 //XXX: Use the psLib function?
    1479 psF32 LinInterpolate(psF32 x0,
    1480                      psF32 x1,
    1481                      psF32 y0,
    1482                      psF32 y1,
    1483                      psF32 y)
    1484 {
    1485     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    1486     psTrace(__func__, 5, "LinInterpolate(%.2f %.2f %.2f %.2f %.2f)\n", x0, x1, y0, y1, y);
    1487     if (y0 == y1) {
    1488         psLogMsg(__func__, PS_LOG_WARN, "WARNING: y0 == y1.  Cannot interpolate.  Returning NaN.\n");
    1489         return(NAN);
    1490     }
    1491 
    1492     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    1493     return(x0 + y * (x1 - x0) / (y1 - y0));
    1494 }
    1495 
    1496 
    1497 
    1498 
    1499 
    1500 
    1501 
    1502 
    1503 
    1504 
    1505 
    1506 
    1507 
    1508 
    1509 
    1510 
    1511 
    1512 
    1513 
    1514 
    15151382
    15161383/******************************************************************************
    … …  
    15321399XXX: Review and ensure that all memory is free'ed at premature exits.
    15331400*****************************************************************************/
    1534 #define INITIAL_NUM_BINS 500.0
     1401#define INITIAL_NUM_BINS 1000.0
    15351402psS32 p_psVectorRobustStats(const psVector* myVector,
    15361403                            const psVector* errors,
    … …  
    15431410        PS_VECTOR_PRINT_F32(myVector);
    15441411    }
    1545     //    psS32 currentId = psMemGetId();
    1546     //    psS32 memLeaks = 0;
    15471412    psHistogram *robustHistogram = NULL;
    15481413    psHistogram *cumulativeRobustHistogram = NULL;
    … …  
    15891454            psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
    15901455            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / INITIAL_NUM_BINS;
    1591             psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
    15921456        }
    15931457        psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize);
    … …  
    16801544        // XXX: Check return codes.
    16811545        //
    1682         printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n");
    16831546        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
    16841547                                  *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    17141577            }
    17151578        }
    1716         psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n", binLo, binHi);
     1579        psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n",
     1580                binLo, binHi);
    17171581        psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo));
    17181582        psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi));
    … …  
    17501614                           binHi,
    17511615                           totalDataPoints * 0.841345f);
    1752             psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
    1753         }
    1754 
     1616            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
     1617                    binLoF32, binHiF32);
     1618        }
     1619
     1620        //
     1621        // This code basically interpolates to find the positions exactly.
     1622        //
    17551623        if (1) {
    1756             psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n", binLo, cumulativeRobustHistogram->nums->data.F32[binLo], cumulativeRobustHistogram->nums->data.F32[binLo+1]);
    1757             psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n", binHi, cumulativeRobustHistogram->nums->data.F32[binHi], cumulativeRobustHistogram->nums->data.F32[binHi+1]);
    1758             // XXX: Add checks to ensure that these bins are not the first (seg fault).
    1759             binLoF32 = LinInterpolate(
    1760                            cumulativeRobustHistogram->bounds->data.F32[binLo],
    1761                            cumulativeRobustHistogram->bounds->data.F32[binLo+1],
    1762                            cumulativeRobustHistogram->nums->data.F32[binLo-1],
    1763                            cumulativeRobustHistogram->nums->data.F32[binLo],
    1764                            totalDataPoints * 0.158655f);
    1765             binHiF32 = LinInterpolate(
    1766                            cumulativeRobustHistogram->bounds->data.F32[binHi],
    1767                            cumulativeRobustHistogram->bounds->data.F32[binHi+1],
    1768                            cumulativeRobustHistogram->nums->data.F32[binHi-1],
    1769                            cumulativeRobustHistogram->nums->data.F32[binHi],
    1770                            totalDataPoints * 0.841345f);
    1771             // XXX: Check for NANs
     1624            psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
     1625                    binLo, cumulativeRobustHistogram->nums->data.F32[binLo],
     1626                    cumulativeRobustHistogram->nums->data.F32[binLo+1]);
     1627            psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
     1628                    binHi, cumulativeRobustHistogram->nums->data.F32[binHi],
     1629                    cumulativeRobustHistogram->nums->data.F32[binHi+1]);
     1630
     1631            psF32 deltaNums;
     1632            psF32 deltaBounds;
     1633            psF32 prevPixels;
     1634            psF32 percentNums;
     1635            psF32 base;
     1636            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binLo+1] - cumulativeRobustHistogram->bounds->data.F32[binLo];
     1637            if (binLo == 0) {
     1638                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
     1639                prevPixels = 0;
     1640            } else {
     1641                deltaNums = cumulativeRobustHistogram->nums->data.F32[binLo] - cumulativeRobustHistogram->nums->data.F32[binLo-1];
     1642                prevPixels = cumulativeRobustHistogram->nums->data.F32[binLo-1];
     1643            }
     1644            percentNums = (totalDataPoints * 0.158655f) - prevPixels;
     1645            base = cumulativeRobustHistogram->bounds->data.F32[binLo];
     1646            binLoF32 = base + (deltaBounds / deltaNums) * percentNums;
     1647            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
     1648                    base, deltaBounds, deltaNums, prevPixels, percentNums);
     1649
     1650            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binHi+1] - cumulativeRobustHistogram->bounds->data.F32[binHi];
     1651            if (binHi == 0) {
     1652                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
     1653                prevPixels = 0;
     1654            } else {
     1655                deltaNums = cumulativeRobustHistogram->nums->data.F32[binHi] - cumulativeRobustHistogram->nums->data.F32[binHi-1];
     1656                prevPixels = cumulativeRobustHistogram->nums->data.F32[binHi-1];
     1657            }
     1658            percentNums = (totalDataPoints * 0.841345f) - prevPixels;
     1659            base = cumulativeRobustHistogram->bounds->data.F32[binHi];
     1660            binHiF32 = base + (deltaBounds / deltaNums) * percentNums;
     1661            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
     1662                    base, deltaBounds, deltaNums, prevPixels, percentNums);
    17721663            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
    17731664        }
    … …  
    17891680        if (sigma < (2 * binSize)) {
    17901681            psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
    1791             psF32 medianLo = robustHistogram->bounds->data.F32[PS_MAX(0, (binMedian - 25))];
    1792             psF32 medianHi = robustHistogram->bounds->data.F32[PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25))];
     1682            psS32 maskLo = PS_MAX(0, (binMedian - 25));
     1683            psS32 maskHi = PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25));
     1684            psF32 medianLo = robustHistogram->bounds->data.F32[maskLo];
     1685            psF32 medianHi = robustHistogram->bounds->data.F32[maskHi];
    17931686            psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
    1794             psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n", binMedian, binMedian - 25, binMedian + 25);
     1687            psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n",
     1688                    binMedian, maskLo, maskHi);
    17951689            psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
    17961690            for (psS32 i = 0 ; i < myVector->n ; i++) {
    17971691                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
    17981692                    tmpMaskVec->data.U8[i] = 1;
    1799                     psTrace(__func__, 8, "Masking element %d is %f\n", i, myVector->data.F32[i]);
     1693                    psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]);
    18001694                }
    18011695            }
    … …  
    18471741    // XXX: Check for errors.
    18481742    //
    1849     printf("ADD: step 8: Interpolate for the lower quartile positions.\n");
    18501743    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
    18511744                           *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    18531746                           binLo25,
    18541747                           totalDataPoints * 0.25f);
    1855     printf("ADD: step 8: Interpolate for the upper quartile positions.\n");
    18561748    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
    18571749                           *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    18911783        psF32 newBinSize = sigma / dN;
    18921784
     1785        //
     1786        // Determine the min/max of the vector (which prior outliers masked out)
     1787        //
    18931788        rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    18941789        rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    … …  
    19031798        }
    19041799
     1800        //
     1801        // Calculate the number of bins.
     1802        //
    19051803        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
    19061804        psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
    … …  
    19861884        }
    19871885        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    1988         psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    19891886        psArray *x = psArrayAlloc((1 + (binMax - binMin)));
    19901887        psS32 j = 0;
    … …  
    19941891            x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
    19951892            ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
    1996             xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);
    19971893            j++;
    19981894        }
    … …  
    20051901        psF32 minY = FLT_MAX;
    20061902        psF32 maxY = -FLT_MAX;
    2007         for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
     1903        for (psS32 i = 0 ; i < y->n ; i++) {
    20081904            if (y->data.F32[i] > maxY) {
    20091905                maxY = y->data.F32[i];
    … …  
    20171913        // XXX: Use the normalize routines for this.
    20181914        //
    2019         for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
     1915        for (psS32 i = 0 ; i < y->n ; i++) {
    20201916            y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
    20211917        }
    … …  
    20311927            PS_VECTOR_PRINT_F32(y);
    20321928        }
    2033         psFree(xTmp);
     1929
     1930        //
     1931        // Fit a Gaussian to the data.
     1932        //
    20341933        rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
    2035 
    20361934        if (rcBool != true) {
    20371935            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    25662464Returns
    25672465    The stats structure.
     2466 
     2467XXX: Should we free stats if the asserts fail?
    25682468 *****************************************************************************/
    25692469psStats* psVectorStats(psStats* stats,
    … …  
    25752475    psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);
    25762476    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    2577     PS_ASSERT_VECTOR_NON_NULL(in, stats);
     2477    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
    25782478    if (mask != NULL) {
    25792479        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);
Note: See TracChangeset for help on using the changeset viewer.