IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 16, 2004, 1:04:18 PM (22 years ago)
Author:
gusciora
Message:

Modified the stdev and mean code.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2724 r2739  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.102 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-16 03:01:40 $
     11 *  @version $Revision: 1.103 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-16 23:04:18 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    11181118 
    11191119XXX: After you fit the polynomial, solve for X analytically.
     1120 
     1121XXX: the vectors do not have to be the same length.  Must insert the proper
     1122tests to ensure that binNum is within acceptable ranges for both vectors.
    11201123*****************************************************************************/
    11211124float fitQuadraticSearchForYThenReturnX(psVector *xVec,
    … …  
    11261129    PS_VECTOR_CHECK_NULL(xVec, NAN);
    11271130    PS_VECTOR_CHECK_NULL(yVec, NAN);
    1128     PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F64, NAN);
    1129     PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F64, NAN);
    1130     PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN);
     1131    PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F32, NAN);
     1132    PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F32, NAN);
     1133    //    PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN);
    11311134    PS_INT_CHECK_RANGE(binNum, 0, (xVec->n - 1), NAN);
    1132 
    1133     PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
    1134     PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
    1135     PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
    1136     PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);
     1135    PS_INT_CHECK_RANGE(binNum, 0, (yVec->n - 1), NAN);
     1136
     1137    //    PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
     1138    //    PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
     1139    //    PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
     1140    //    PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);
     1141    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1142    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1143    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1144    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
     1145
    11371146    float tmpFloat;
    11381147
    … …  
    11621171                    false,
    11631172                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
     1173            psFree(myPoly);
     1174            psFree(x);
     1175            psFree(y);
     1176            psFree(yErr);
    11641177            return(NAN);
    11651178        }
    … …  
    11711184                    false,
    11721185                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
     1186            psFree(myPoly);
     1187            psFree(x);
     1188            psFree(y);
     1189            psFree(yErr);
    11731190            return(NAN);
    11741191        }
    … …  
    11921209    }
    11931210
     1211    psFree(myPoly);
     1212    psFree(x);
     1213    psFree(y);
     1214    psFree(yErr);
    11941215    return(tmpFloat);
    11951216}
    … …  
    13941415    myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    13951416
    1396 
    13971417    // Using the above (myMean, myStdev) as initial estimates, we fit a
    13981418    // Gaussian to the robustHistogramVector.
    … …  
    14111431    myParams->data.F32[0] = myMean;
    14121432    myParams->data.F32[1] = myStdev;
     1433    psImage *covar = NULL;
    14131434    psMinimizeLMChi2(min,
    14141435                     NULL,
    … …  
    14191440                     NULL,
    14201441                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
     1442    psFree(covar);
    14211443    psFree(min);
    14221444    psFree(myParams);
    … …  
    14291451        if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
    14301452            psLogMsg(__func__, PS_LOG_WARN,
    1431                      "WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
     1453                     "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
    14321454            psLogMsg(__func__, PS_LOG_WARN,
    1433                      "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n",
     1455                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
    14341456                     myMean, myParams->data.F32[0]);
    14351457            stats->robustMean = myMean;
    … …  
    14461468        if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
    14471469            psLogMsg(__func__, PS_LOG_WARN,
    1448                      "WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
     1470                     "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
    14491471            psLogMsg(__func__, PS_LOG_WARN,
    1450                      "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n",
    1451                      myMean, myParams->data.F32[0]);
     1472                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
     1473                     myStdev, myParams->data.F32[1]);
    14521474            stats->robustStdev = myStdev;
    14531475        } else {
    … …  
    14631485        // Take a psVector.  Fit a polynomial y = f(x) to the 3 data elements
    14641486        // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0.
    1465         stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
    1466                               cumulativeRobustSums,
    1467                               medianBinNum,
    1468                               sumRobust/2.0);
     1487
     1488        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
     1489        //Determine id that is okay.
     1490        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
     1491                                  robustHistogram->bounds,
     1492                                  cumulativeRobustSums,
     1493                                  medianBinNum,
     1494                                  sumRobust/2.0);
    14691495    }
    14701496
    … …  
    14761502        countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
    14771503
    1478         stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
    1479                           cumulativeRobustSums,
    1480                           LQBinNum,
    1481                           countFloat/4.0);
    1482         stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
    1483                           cumulativeRobustSums,
    1484                           UQBinNum,
    1485                           3.0 * countFloat/4.0);
     1504        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
     1505        //Determine id that is okay.
     1506        stats->robustLQ = fitQuadraticSearchForYThenReturnX(
     1507                              robustHistogram->bounds,
     1508                              cumulativeRobustSums,
     1509                              LQBinNum,
     1510                              countFloat/4.0);
     1511        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
     1512        //Determine id that is okay.
     1513        stats->robustUQ = fitQuadraticSearchForYThenReturnX(
     1514                              robustHistogram->bounds,
     1515                              cumulativeRobustSums,
     1516                              UQBinNum,
     1517                              3.0 * countFloat/4.0);
    14861518    }
    14871519    // XXX: I think sumNfit == sumN50 here.
Note: See TracChangeset for help on using the changeset viewer.