IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2269


Ignore:
Timestamp:
Nov 2, 2004, 5:30:30 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFunctions.h

    r2204 r2269  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-10-27 00:57:31 $
     14*  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-11-03 03:30:30 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                               const psSpline1D *spline);
    410410
    411 psS32 p_psVectorBinDisectF32(float *bins,
    412                              psS32 numBins,
    413                              float x);
    414 
    415 psS32 p_psVectorBinDisectS32(psS32 *bins,
    416                              psS32 numBins,
    417                              psS32 x);
    418 
    419411psS32 p_psVectorBinDisect(psVector *bins,
    420412                          psScalar *x);

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

    r2267 r2269  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-11-02 19:08:33 $
     11 *  @version $Revision: 1.84 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-11-03 03:30:30 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6161returned as a psVector sums.
    6262 
    63 XXX: change name
    6463XXX: Use a static vector.
    6564 *****************************************************************************/
    66 void p_psBuildSums1D(double x,
    67                      psS32 polyOrder,
    68                      psVector* sums)
     65void psBuildSums1D(double x,
     66                   psS32 polyOrder,
     67                   psVector* sums)
    6968{
    7069    psS32 i = 0;
    … …  
    7473        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
    7574    }
    76 
    7775    if (polyOrder > sums->n) {
    7876        sums = psVectorRealloc(sums, polyOrder);
    … …  
    109107    psS32 i;
    110108    psS32 k;
    111     psBool mustFreeX = false;
    112109    float sig;
    113110    float p;
    … …  
    150147    }
    151148
    152     if (mustFreeX == true) {
    153         psFree(X);
    154     }
    155149    psFree(u);
    156 
    157150    psTrace(".psLib.dataManip.CalculateSecondDerivs", 4,
    158151            "---- CalculateSecondDerivs() end ----\n");
    … …  
    171164polynomials which are stored in psSpline1D.
    172165 
    173 XXX: check types/sizes
    174166XXX: This is F32 only
    175167 *****************************************************************************/
     168/*
    176169float p_psNRSpline1DEval(psSpline1D *spline,
    177170                         const psVector* restrict x,
    … …  
    187180    PS_VECTOR_CHECK_NULL(y, NAN);
    188181    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NAN);
    189 
     182 
    190183    psS32 n;
    191184    psS32 klo;
    … …  
    197190    float D;
    198191    float Y;
    199 
     192 
    200193    n = spline->n;
    201     klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X);
     194    klo = p_psVectorBinDisect32(spline->domains, (spline->n)+1, X);
    202195    khi = klo + 1;
    203196    H = (spline->domains)[khi] - (spline->domains)[klo];
    … …  
    206199    C = ((A*A*A)-A) * (H*H/6.0);
    207200    D = ((B*B*B)-B) * (H*H/6.0);
    208 
     201 
    209202    Y = (A * y->data.F32[klo]) +
    210203        (B * y->data.F32[khi]) +
    211204        (C * (spline->p_psDeriv2)[klo]) +
    212205        (D * (spline->p_psDeriv2)[khi]);
    213 
     206 
    214207    return(Y);
    215208}
    216 
     209*/
    217210/*****************************************************************************/
    218211/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    239232 
    240233XXX: usage of yErr is not specified in IfA documentation.
     234 
    241235XXX: Is the x argument redundant?  What do we do if the x argument is
    242236supplied, but does not equal the domains specified in mySpline?
    … …  
    414408/******************************************************************************
    415409XXX: We assume unnormalized gaussians.
    416 XXX: Currently, yErr is ignored.
    417410 *****************************************************************************/
    418411psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
    … …  
    441434        x = ((psVector *) (coords->data[i]))->data.F32[0];
    442435        out->data.F32[i] = psGaussian(x, mean, stdev, false);
    443         //        psTrace(".psLib.dataManip.psMinimize", 6,
    444         //            "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]);
    445436    }
    446437
    … …  
    449440        float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    450441        deriv->data.F32[i][0] = tmp / (stdev * stdev);
    451 
    452442        tmp = (x - mean) * (x - mean) *
    453443              psGaussian(x, mean, stdev, 0);
    … …  
    482472XXX: This must work for both F32 and F64.  F32 is currently implemented.
    483473     Note: since the LUD routines are only implemented in F64, then we
    484 >     will have to convert all F32 input vectors to F64 regardless.  So,
     474     will have to convert all F32 input vectors to F64 regardless.  So,
    485475     the F64 port might be.
    486476 *****************************************************************************/
    … …  
    730720
    731721/******************************************************************************
    732 p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
    733 polynomial of degree myPoly to the data points (x, y) and return the
    734 coefficients of that polynomial.
     722VectorFitPolynomial1DCheb():  This routine will fit a Chebyshev polynomial of
     723degree myPoly to the data points (x, y) and return the coefficients of that
     724polynomial.
    735725 
    736726XXX: yErr is currently ignored.
    … …  
    738728XXX: Use private name?
    739729*****************************************************************************/
    740 psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
     730psPolynomial1D *VectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
    741731        const psVector* restrict x,
    742732        const psVector* restrict y,
    … …  
    785775        psFree(fScalar);
    786776
    787         psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6,
     777        psTrace(".psLib.dataManip.VectorFitPolynomial1DCheby", 6,
    788778                "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
    789779                x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
    … …  
    809799
    810800/******************************************************************************
    811 p_psVectorFitPolynomial1DOrd():  This routine will fit an ordinary
    812 polynomial of degree myPoly to the data points (x, y) and return the
    813 coefficients of that polynomial.
     801VectorFitPolynomial1DOrd():  This routine will fit an ordinary polynomial of
     802degree myPoly to the data points (x, y) and return the coefficients of that
     803polynomial.
    814804 
    815805XXX: Use private name?
    816806XXX: Use recycled vectors.
    817807 *****************************************************************************/
    818 psPolynomial1D* p_psVectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
     808psPolynomial1D* VectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
    819809        const psVector* restrict x,
    820810        const psVector* restrict y,
    … …  
    833823    psVector* xSums = NULL;
    834824
    835     psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
    836             "---- psVectorFitPolynomial1D() begin ----\n");
    837     // printf("psVectorFitPolynomial1D()\n");
     825    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
     826            "---- VectorFitPolynomial1DOrd() begin ----\n");
     827    // printf("VectorFitPolynomial1D()\n");
    838828    // for (i=0;i<x->n;i++) {
    839829    // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
    … …  
    865855    if (yErr == NULL) {
    866856        for (i = 0; i < X->n; i++) {
    867             p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
     857            psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
    868858
    869859            for (k = 0; k < polyOrder; k++) {
    … …  
    879869    } else {
    880870        for (i = 0; i < X->n; i++) {
    881             p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
     871            psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
    882872
    883873            for (k = 0; k < polyOrder; k++) {
    … …  
    916906    psFree(xSums);
    917907
    918     psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
    919             "---- psVectorFitPolynomial1D() begin ----\n");
     908    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
     909            "---- VectorFitPolynomial1DOrd() begin ----\n");
    920910    return (myPoly);
    921911}
    … …  
    962952        PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n);
    963953        if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    964             p_psNormalizeVector(x64Static);
     954            p_psNormalizeVectorRange(x64Static, -1.0, 1.0);
    965955        }
    966956        x64 = x64Static;
    … …  
    974964    // Call the appropriate vector fitting routine.
    975965    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    976         p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     966        VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    977967    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    978         tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
     968        tmpPoly = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    979969    } else {
    980970        psError(__func__, "unknown polynomial type.\n");

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

    r2268 r2269  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-11-02 19:13:03 $
     11 *  @version $Revision: 1.84 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-11-03 03:30:30 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    403403 *****************************************************************************/
    404404void p_psVectorSampleMedian(const psVector* restrict myVector,
    405                             const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     405                            const psVector* restrict maskVector,
     406                            psU32 maskVal,
     407                            psStats* stats)
    406408{
    407409    psVector* unsortedVector = NULL;    // Temporary vector
    408410    psVector* sortedVector = NULL;      // Temporary vector
    409     psS32 i = 0;                  // Loop index variable
    410     psS32 count = 0;              // # of points in this mean?
    411     psS32 nValues = 0;            // # of points in vector
    412     float rangeMin = 0.0;       // Exclude data below this
    413     float rangeMax = 0.0;       // Exclude date above this
     411    psS32 i = 0;                        // Loop index variable
     412    psS32 count = 0;                    // # of points in this mean?
     413    psS32 nValues = 0;                  // # of points in vector
     414    float rangeMin = 0.0;               // Exclude data below this
     415    float rangeMax = 0.0;               // Exclude date above this
    414416
    415417    // Determine how many data points fit inside this min/max range
    416     // and are not masked, IF the maskVector is not NULL>
     418    // and are not masked, if the maskVector is not NULL>
    417419    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    418420
    … …  
    551553 *****************************************************************************/
    552554void p_psVectorSampleQuartiles(const psVector* restrict myVector,
    553                                const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     555                               const psVector* restrict maskVector,
     556                               psU32 maskVal,
     557                               psStats* stats)
    554558{
    555559    psVector* unsortedVector = NULL;    // Temporary vector
    556560    psVector* sortedVector = NULL;      // Temporary vector
    557561    psS32 i = 0;                  // Loop index variable
    558     psS32 count = 0;              // # of points in this mean?
     562    psS32 count = 0;              // # of points in this mean.
    559563    psS32 nValues = 0;            // # data points
    560564    float rangeMin = 0.0;       // Exclude data below this
    … …  
    631635 *****************************************************************************/
    632636void p_psVectorSampleStdev(const psVector* restrict myVector,
    633                            const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     637                           const psVector* restrict maskVector,
     638                           psU32 maskVal,
     639                           psStats* stats)
    634640{
    635641    psS32 i = 0;                  // Loop index variable
    … …  
    713719    stats
    714720Returns
    715     NULL
    716  *****************************************************************************/
    717 void p_psVectorClippedStats(const psVector* restrict myVector,
    718                             const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     721    0 for success.
     722 *****************************************************************************/
     723int p_psVectorClippedStats(const psVector* restrict myVector,
     724                           const psVector* restrict maskVector,
     725                           psU32 maskVal,
     726                           psStats* stats)
    719727{
    720728    psS32 i = 0;                  // Loop index variable
    … …  
    726734    psVector* tmpMask = NULL;   // Temporary vector
    727735
    728     // Endure that stats->clipIter is within the proper range.
     736    // Ensure that stats->clipIter is within the proper range.
    729737    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
    730738        psError(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
    731     }
    732     // Endure that stats->clipSigma is within the proper range.
     739        return(-1);
     740    }
     741    // Ensure that stats->clipSigma is within the proper range.
    733742    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
    734743        psError(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
     744        return(-1);
    735745    }
    736746    // We allocate a temporary mask vector since during the iterative
    … …  
    798808
    799809    psFree(tmpMask);
    800 }
    801 
    802 void p_psNormalizeVectorF32_0(psVector* myData)
     810    return(0);
     811}
     812
     813/*****************************************************************************
     814 *****************************************************************************/
     815void p_psNormalizeVectorRangeF32(psVector* myData,
     816                                 float outLow,
     817                                 float outHigh)
    803818{
    804819    float min = (float)HUGE;
    … …  
    815830    }
    816831
    817     //  float range = max - min;
    818     //    for (i = 0; i < myData->n; i++) {
    819     //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
    820     //    }
    821832    for (i = 0; i < myData->n; i++) {
    822         myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min);
    823     }
    824 }
    825 
    826 
    827 
    828 /*****************************************************************************
    829 p_psNormalizeVectorF32(myData): this is a private function which normalizes the
    830 elements of a vector to a range between 0.0 and 1.0.
    831  
    832 XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)?
    833  *****************************************************************************/
    834 void p_psNormalizeVectorF32(psVector* myData)
     833        myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
     834    }
     835}
     836void p_psNormalizeVectorRangeF64(psVector* myData,
     837                                 float outLow,
     838                                 float outHigh)
    835839{
    836840    float min = (float)HUGE;
    … …  
    840844    for (i = 0; i < myData->n; i++) {
    841845        if (myData->data.F32[i] < min) {
    842             min = myData->data.F32[i];
     846            min = myData->data.F64[i];
    843847        }
    844848        if (myData->data.F32[i] > max) {
    845             max = myData->data.F32[i];
    846         }
    847     }
    848 
    849     //  float range = max - min;
    850     //    for (i = 0; i < myData->n; i++) {
    851     //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
    852     //    }
     849            max = myData->data.F64[i];
     850        }
     851    }
     852
    853853    for (i = 0; i < myData->n; i++) {
    854         myData->data.F32[i] = (myData->data.F32[i] - 0.5 * (min + max)) /
    855                               (0.5 * (max - min));
     854        myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
    856855    }
    857856}
    858857
    859858/*****************************************************************************
    860 p_psNormalizeVectorF64(myData): this is a private function which normalizes the
    861 elements of a vector to a range between -1.0 and 1.0.
    862 XXX: 0-1 or -1:1?
    863  
    864 XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)?
    865  *****************************************************************************/
    866 void p_psNormalizeVectorF64(psVector* myData)
    867 {
    868     float min = (double)HUGE;
    869     float max = (double)-HUGE;
    870     double range = 0.0;
    871     psS32 i = 0;
    872 
    873     for (i = 0; i < myData->n; i++) {
    874         if (myData->data.F64[i] < min) {
    875             min = myData->data.F64[i];
    876         }
    877         if (myData->data.F64[i] > max) {
    878             max = myData->data.F64[i];
    879         }
    880     }
    881 
    882     range = max - min;
    883     for (i = 0; i < myData->n; i++) {
    884         myData->data.F64[i] = -1.0 + 2.0 *
    885                               ((myData->data.F64[i] - min) / range);
    886     }
    887 
    888     //    for (i = 0; i < myData->n; i++) {
    889     //        myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) /
    890     //                              (0.5 * (max - min));
    891     //    }
    892 }
    893 
    894 // XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)?
    895 void p_psNormalizeVector(psVector* myData)
     859psNormalizeVectorRange(myData, low, high): this is a private function which
     860normalizes the elements of a vector to a range between low and high.
     861 *****************************************************************************/
     862void psNormalizeVectorRange(psVector* myData,
     863                            float low,
     864                            float high)
    896865{
    897866    if (myData->type.type == PS_TYPE_F32) {
    898         p_psNormalizeVectorF32(myData);
     867        p_psNormalizeVectorRangeF32(myData, low, high);
    899868    } else if (myData->type.type == PS_TYPE_F64) {
    900         p_psNormalizeVectorF64(myData);
     869        p_psNormalizeVectorRangeF64(myData, low, high);
    901870    } else {
    902871        psError(__func__, "Unalowable data type.\n");
    … …  
    914883 
    915884XXX: Terminate when f(x)-getThisValue is within some error tolerance.
     885 
     886XXX: Solve for X analytically.
    916887 *****************************************************************************/
    917888float p_ps1DPolyMedian(psPolynomial1D* myPoly,
    … …  
    924895    float oldMidpoint = 1.0;
    925896    float f = 0.0;
    926 
    927     // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
    928897
    929898    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
    … …  
    956925and i+1 is used for x[i]).  It then determines for what value x does that
    957926quadratic f(x) = yVal (the input parameter).
     927 
     928XXX: After you fit the polynomial, solve for X analytically.
    958929*****************************************************************************/
    959930float fitQuadraticSearchForYThenReturnX(psVector *xVec,
    … …  
    10341005    stats
    10351006Returns
    1036     NULL
     1007    0 on success.
    10371008*****************************************************************************/
    10381009int p_psVectorRobustStats(const psVector* restrict myVector,
    … …  
    11191090                            tmpStats->clippedStdev / 4.0f);
    11201091
    1121     // The following was necessary to fit a gaussian to the data, since
    1122     // gaussian functions produce data between 0.0 and 1.0.
    1123     // p_psNormalizeVectorF32(robustHistogramVector);
    1124 
    11251092    /**************************************************************************
    11261093    Determine the median/lower/upper quartile bin numbers.
    … …  
    12071174    psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    12081175
    1209     p_psNormalizeVectorF32_0(robustHistogramVector);
     1176    p_psNormalizeVectorRange(robustHistogramVector, 0.0, 1.0);
    12101177    for (i=0;i<robustHistogramVector->n;i++) {
    12111178        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    … …  
    16201587    }
    16211588    // ************************************************************************
    1622     // XXX: The Stdev calculation requires the mean.  Should we assume the
    1623     // mean has already been calculated? Or should we always calculate it?
    16241589    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    16251590        p_psVectorSampleMean(inF32, mask, maskVal, stats);
    … …  
    16431608
    16441609    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    1645         p_psVectorClippedStats(inF32, mask, maskVal, stats);
     1610        if (0 != p_psVectorClippedStats(inF32, mask, maskVal, stats)) {
     1611            psError(__func__, "p_psVectorClippedStats() failed.\n");
     1612        }
    16461613    }
    16471614    // ************************************************************************

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

    r2204 r2269  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-10-27 00:57:31 $
     12 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-11-03 03:30:30 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    174174);
    175175
    176 
    177 void p_psNormalizeVector(psVector* myData);
    178 
    179 void p_psNormalizeVectorF64(psVector* myData);
     176void psNormalizeVectorRange(psVector* myData,
     177                            float low,
     178                            float high);
    180179
    181180/// @}

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

    r2267 r2269  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-11-02 19:08:33 $
     11 *  @version $Revision: 1.84 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-11-03 03:30:30 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6161returned as a psVector sums.
    6262 
    63 XXX: change name
    6463XXX: Use a static vector.
    6564 *****************************************************************************/
    66 void p_psBuildSums1D(double x,
    67                      psS32 polyOrder,
    68                      psVector* sums)
     65void psBuildSums1D(double x,
     66                   psS32 polyOrder,
     67                   psVector* sums)
    6968{
    7069    psS32 i = 0;
    … …  
    7473        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
    7574    }
    76 
    7775    if (polyOrder > sums->n) {
    7876        sums = psVectorRealloc(sums, polyOrder);
    … …  
    109107    psS32 i;
    110108    psS32 k;
    111     psBool mustFreeX = false;
    112109    float sig;
    113110    float p;
    … …  
    150147    }
    151148
    152     if (mustFreeX == true) {
    153         psFree(X);
    154     }
    155149    psFree(u);
    156 
    157150    psTrace(".psLib.dataManip.CalculateSecondDerivs", 4,
    158151            "---- CalculateSecondDerivs() end ----\n");
    … …  
    171164polynomials which are stored in psSpline1D.
    172165 
    173 XXX: check types/sizes
    174166XXX: This is F32 only
    175167 *****************************************************************************/
     168/*
    176169float p_psNRSpline1DEval(psSpline1D *spline,
    177170                         const psVector* restrict x,
    … …  
    187180    PS_VECTOR_CHECK_NULL(y, NAN);
    188181    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NAN);
    189 
     182 
    190183    psS32 n;
    191184    psS32 klo;
    … …  
    197190    float D;
    198191    float Y;
    199 
     192 
    200193    n = spline->n;
    201     klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X);
     194    klo = p_psVectorBinDisect32(spline->domains, (spline->n)+1, X);
    202195    khi = klo + 1;
    203196    H = (spline->domains)[khi] - (spline->domains)[klo];
    … …  
    206199    C = ((A*A*A)-A) * (H*H/6.0);
    207200    D = ((B*B*B)-B) * (H*H/6.0);
    208 
     201 
    209202    Y = (A * y->data.F32[klo]) +
    210203        (B * y->data.F32[khi]) +
    211204        (C * (spline->p_psDeriv2)[klo]) +
    212205        (D * (spline->p_psDeriv2)[khi]);
    213 
     206 
    214207    return(Y);
    215208}
    216 
     209*/
    217210/*****************************************************************************/
    218211/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    239232 
    240233XXX: usage of yErr is not specified in IfA documentation.
     234 
    241235XXX: Is the x argument redundant?  What do we do if the x argument is
    242236supplied, but does not equal the domains specified in mySpline?
    … …  
    414408/******************************************************************************
    415409XXX: We assume unnormalized gaussians.
    416 XXX: Currently, yErr is ignored.
    417410 *****************************************************************************/
    418411psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
    … …  
    441434        x = ((psVector *) (coords->data[i]))->data.F32[0];
    442435        out->data.F32[i] = psGaussian(x, mean, stdev, false);
    443         //        psTrace(".psLib.dataManip.psMinimize", 6,
    444         //            "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]);
    445436    }
    446437
    … …  
    449440        float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    450441        deriv->data.F32[i][0] = tmp / (stdev * stdev);
    451 
    452442        tmp = (x - mean) * (x - mean) *
    453443              psGaussian(x, mean, stdev, 0);
    … …  
    482472XXX: This must work for both F32 and F64.  F32 is currently implemented.
    483473     Note: since the LUD routines are only implemented in F64, then we
    484 >     will have to convert all F32 input vectors to F64 regardless.  So,
     474     will have to convert all F32 input vectors to F64 regardless.  So,
    485475     the F64 port might be.
    486476 *****************************************************************************/
    … …  
    730720
    731721/******************************************************************************
    732 p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
    733 polynomial of degree myPoly to the data points (x, y) and return the
    734 coefficients of that polynomial.
     722VectorFitPolynomial1DCheb():  This routine will fit a Chebyshev polynomial of
     723degree myPoly to the data points (x, y) and return the coefficients of that
     724polynomial.
    735725 
    736726XXX: yErr is currently ignored.
    … …  
    738728XXX: Use private name?
    739729*****************************************************************************/
    740 psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
     730psPolynomial1D *VectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
    741731        const psVector* restrict x,
    742732        const psVector* restrict y,
    … …  
    785775        psFree(fScalar);
    786776
    787         psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6,
     777        psTrace(".psLib.dataManip.VectorFitPolynomial1DCheby", 6,
    788778                "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
    789779                x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
    … …  
    809799
    810800/******************************************************************************
    811 p_psVectorFitPolynomial1DOrd():  This routine will fit an ordinary
    812 polynomial of degree myPoly to the data points (x, y) and return the
    813 coefficients of that polynomial.
     801VectorFitPolynomial1DOrd():  This routine will fit an ordinary polynomial of
     802degree myPoly to the data points (x, y) and return the coefficients of that
     803polynomial.
    814804 
    815805XXX: Use private name?
    816806XXX: Use recycled vectors.
    817807 *****************************************************************************/
    818 psPolynomial1D* p_psVectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
     808psPolynomial1D* VectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
    819809        const psVector* restrict x,
    820810        const psVector* restrict y,
    … …  
    833823    psVector* xSums = NULL;
    834824
    835     psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
    836             "---- psVectorFitPolynomial1D() begin ----\n");
    837     // printf("psVectorFitPolynomial1D()\n");
     825    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
     826            "---- VectorFitPolynomial1DOrd() begin ----\n");
     827    // printf("VectorFitPolynomial1D()\n");
    838828    // for (i=0;i<x->n;i++) {
    839829    // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
    … …  
    865855    if (yErr == NULL) {
    866856        for (i = 0; i < X->n; i++) {
    867             p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
     857            psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
    868858
    869859            for (k = 0; k < polyOrder; k++) {
    … …  
    879869    } else {
    880870        for (i = 0; i < X->n; i++) {
    881             p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
     871            psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
    882872
    883873            for (k = 0; k < polyOrder; k++) {
    … …  
    916906    psFree(xSums);
    917907
    918     psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
    919             "---- psVectorFitPolynomial1D() begin ----\n");
     908    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
     909            "---- VectorFitPolynomial1DOrd() begin ----\n");
    920910    return (myPoly);
    921911}
    … …  
    962952        PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n);
    963953        if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    964             p_psNormalizeVector(x64Static);
     954            p_psNormalizeVectorRange(x64Static, -1.0, 1.0);
    965955        }
    966956        x64 = x64Static;
    … …  
    974964    // Call the appropriate vector fitting routine.
    975965    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    976         p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     966        VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    977967    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    978         tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
     968        tmpPoly = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    979969    } else {
    980970        psError(__func__, "unknown polynomial type.\n");

  • trunk/psLib/src/math/psPolynomial.h

    r2204 r2269  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-10-27 00:57:31 $
     14*  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-11-03 03:30:30 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                               const psSpline1D *spline);
    410410
    411 psS32 p_psVectorBinDisectF32(float *bins,
    412                              psS32 numBins,
    413                              float x);
    414 
    415 psS32 p_psVectorBinDisectS32(psS32 *bins,
    416                              psS32 numBins,
    417                              psS32 x);
    418 
    419411psS32 p_psVectorBinDisect(psVector *bins,
    420412                          psScalar *x);

  • trunk/psLib/src/math/psSpline.h

    r2204 r2269  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-10-27 00:57:31 $
     14*  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-11-03 03:30:30 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    409409                               const psSpline1D *spline);
    410410
    411 psS32 p_psVectorBinDisectF32(float *bins,
    412                              psS32 numBins,
    413                              float x);
    414 
    415 psS32 p_psVectorBinDisectS32(psS32 *bins,
    416                              psS32 numBins,
    417                              psS32 x);
    418 
    419411psS32 p_psVectorBinDisect(psVector *bins,
    420412                          psScalar *x);

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

    r2268 r2269  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-11-02 19:13:03 $
     11 *  @version $Revision: 1.84 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-11-03 03:30:30 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    403403 *****************************************************************************/
    404404void p_psVectorSampleMedian(const psVector* restrict myVector,
    405                             const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     405                            const psVector* restrict maskVector,
     406                            psU32 maskVal,
     407                            psStats* stats)
    406408{
    407409    psVector* unsortedVector = NULL;    // Temporary vector
    408410    psVector* sortedVector = NULL;      // Temporary vector
    409     psS32 i = 0;                  // Loop index variable
    410     psS32 count = 0;              // # of points in this mean?
    411     psS32 nValues = 0;            // # of points in vector
    412     float rangeMin = 0.0;       // Exclude data below this
    413     float rangeMax = 0.0;       // Exclude date above this
     411    psS32 i = 0;                        // Loop index variable
     412    psS32 count = 0;                    // # of points in this mean?
     413    psS32 nValues = 0;                  // # of points in vector
     414    float rangeMin = 0.0;               // Exclude data below this
     415    float rangeMax = 0.0;               // Exclude date above this
    414416
    415417    // Determine how many data points fit inside this min/max range
    416     // and are not masked, IF the maskVector is not NULL>
     418    // and are not masked, if the maskVector is not NULL>
    417419    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
    418420
    … …  
    551553 *****************************************************************************/
    552554void p_psVectorSampleQuartiles(const psVector* restrict myVector,
    553                                const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     555                               const psVector* restrict maskVector,
     556                               psU32 maskVal,
     557                               psStats* stats)
    554558{
    555559    psVector* unsortedVector = NULL;    // Temporary vector
    556560    psVector* sortedVector = NULL;      // Temporary vector
    557561    psS32 i = 0;                  // Loop index variable
    558     psS32 count = 0;              // # of points in this mean?
     562    psS32 count = 0;              // # of points in this mean.
    559563    psS32 nValues = 0;            // # data points
    560564    float rangeMin = 0.0;       // Exclude data below this
    … …  
    631635 *****************************************************************************/
    632636void p_psVectorSampleStdev(const psVector* restrict myVector,
    633                            const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     637                           const psVector* restrict maskVector,
     638                           psU32 maskVal,
     639                           psStats* stats)
    634640{
    635641    psS32 i = 0;                  // Loop index variable
    … …  
    713719    stats
    714720Returns
    715     NULL
    716  *****************************************************************************/
    717 void p_psVectorClippedStats(const psVector* restrict myVector,
    718                             const psVector* restrict maskVector, psU32 maskVal, psStats* stats)
     721    0 for success.
     722 *****************************************************************************/
     723int p_psVectorClippedStats(const psVector* restrict myVector,
     724                           const psVector* restrict maskVector,
     725                           psU32 maskVal,
     726                           psStats* stats)
    719727{
    720728    psS32 i = 0;                  // Loop index variable
    … …  
    726734    psVector* tmpMask = NULL;   // Temporary vector
    727735
    728     // Endure that stats->clipIter is within the proper range.
     736    // Ensure that stats->clipIter is within the proper range.
    729737    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
    730738        psError(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
    731     }
    732     // Endure that stats->clipSigma is within the proper range.
     739        return(-1);
     740    }
     741    // Ensure that stats->clipSigma is within the proper range.
    733742    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
    734743        psError(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
     744        return(-1);
    735745    }
    736746    // We allocate a temporary mask vector since during the iterative
    … …  
    798808
    799809    psFree(tmpMask);
    800 }
    801 
    802 void p_psNormalizeVectorF32_0(psVector* myData)
     810    return(0);
     811}
     812
     813/*****************************************************************************
     814 *****************************************************************************/
     815void p_psNormalizeVectorRangeF32(psVector* myData,
     816                                 float outLow,
     817                                 float outHigh)
    803818{
    804819    float min = (float)HUGE;
    … …  
    815830    }
    816831
    817     //  float range = max - min;
    818     //    for (i = 0; i < myData->n; i++) {
    819     //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
    820     //    }
    821832    for (i = 0; i < myData->n; i++) {
    822         myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min);
    823     }
    824 }
    825 
    826 
    827 
    828 /*****************************************************************************
    829 p_psNormalizeVectorF32(myData): this is a private function which normalizes the
    830 elements of a vector to a range between 0.0 and 1.0.
    831  
    832 XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)?
    833  *****************************************************************************/
    834 void p_psNormalizeVectorF32(psVector* myData)
     833        myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
     834    }
     835}
     836void p_psNormalizeVectorRangeF64(psVector* myData,
     837                                 float outLow,
     838                                 float outHigh)
    835839{
    836840    float min = (float)HUGE;
    … …  
    840844    for (i = 0; i < myData->n; i++) {
    841845        if (myData->data.F32[i] < min) {
    842             min = myData->data.F32[i];
     846            min = myData->data.F64[i];
    843847        }
    844848        if (myData->data.F32[i] > max) {
    845             max = myData->data.F32[i];
    846         }
    847     }
    848 
    849     //  float range = max - min;
    850     //    for (i = 0; i < myData->n; i++) {
    851     //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
    852     //    }
     849            max = myData->data.F64[i];
     850        }
     851    }
     852
    853853    for (i = 0; i < myData->n; i++) {
    854         myData->data.F32[i] = (myData->data.F32[i] - 0.5 * (min + max)) /
    855                               (0.5 * (max - min));
     854        myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
    856855    }
    857856}
    858857
    859858/*****************************************************************************
    860 p_psNormalizeVectorF64(myData): this is a private function which normalizes the
    861 elements of a vector to a range between -1.0 and 1.0.
    862 XXX: 0-1 or -1:1?
    863  
    864 XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)?
    865  *****************************************************************************/
    866 void p_psNormalizeVectorF64(psVector* myData)
    867 {
    868     float min = (double)HUGE;
    869     float max = (double)-HUGE;
    870     double range = 0.0;
    871     psS32 i = 0;
    872 
    873     for (i = 0; i < myData->n; i++) {
    874         if (myData->data.F64[i] < min) {
    875             min = myData->data.F64[i];
    876         }
    877         if (myData->data.F64[i] > max) {
    878             max = myData->data.F64[i];
    879         }
    880     }
    881 
    882     range = max - min;
    883     for (i = 0; i < myData->n; i++) {
    884         myData->data.F64[i] = -1.0 + 2.0 *
    885                               ((myData->data.F64[i] - min) / range);
    886     }
    887 
    888     //    for (i = 0; i < myData->n; i++) {
    889     //        myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) /
    890     //                              (0.5 * (max - min));
    891     //    }
    892 }
    893 
    894 // XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)?
    895 void p_psNormalizeVector(psVector* myData)
     859psNormalizeVectorRange(myData, low, high): this is a private function which
     860normalizes the elements of a vector to a range between low and high.
     861 *****************************************************************************/
     862void psNormalizeVectorRange(psVector* myData,
     863                            float low,
     864                            float high)
    896865{
    897866    if (myData->type.type == PS_TYPE_F32) {
    898         p_psNormalizeVectorF32(myData);
     867        p_psNormalizeVectorRangeF32(myData, low, high);
    899868    } else if (myData->type.type == PS_TYPE_F64) {
    900         p_psNormalizeVectorF64(myData);
     869        p_psNormalizeVectorRangeF64(myData, low, high);
    901870    } else {
    902871        psError(__func__, "Unalowable data type.\n");
    … …  
    914883 
    915884XXX: Terminate when f(x)-getThisValue is within some error tolerance.
     885 
     886XXX: Solve for X analytically.
    916887 *****************************************************************************/
    917888float p_ps1DPolyMedian(psPolynomial1D* myPoly,
    … …  
    924895    float oldMidpoint = 1.0;
    925896    float f = 0.0;
    926 
    927     // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
    928897
    929898    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
    … …  
    956925and i+1 is used for x[i]).  It then determines for what value x does that
    957926quadratic f(x) = yVal (the input parameter).
     927 
     928XXX: After you fit the polynomial, solve for X analytically.
    958929*****************************************************************************/
    959930float fitQuadraticSearchForYThenReturnX(psVector *xVec,
    … …  
    10341005    stats
    10351006Returns
    1036     NULL
     1007    0 on success.
    10371008*****************************************************************************/
    10381009int p_psVectorRobustStats(const psVector* restrict myVector,
    … …  
    11191090                            tmpStats->clippedStdev / 4.0f);
    11201091
    1121     // The following was necessary to fit a gaussian to the data, since
    1122     // gaussian functions produce data between 0.0 and 1.0.
    1123     // p_psNormalizeVectorF32(robustHistogramVector);
    1124 
    11251092    /**************************************************************************
    11261093    Determine the median/lower/upper quartile bin numbers.
    … …  
    12071174    psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    12081175
    1209     p_psNormalizeVectorF32_0(robustHistogramVector);
     1176    p_psNormalizeVectorRange(robustHistogramVector, 0.0, 1.0);
    12101177    for (i=0;i<robustHistogramVector->n;i++) {
    12111178        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    … …  
    16201587    }
    16211588    // ************************************************************************
    1622     // XXX: The Stdev calculation requires the mean.  Should we assume the
    1623     // mean has already been calculated? Or should we always calculate it?
    16241589    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    16251590        p_psVectorSampleMean(inF32, mask, maskVal, stats);
    … …  
    16431608
    16441609    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    1645         p_psVectorClippedStats(inF32, mask, maskVal, stats);
     1610        if (0 != p_psVectorClippedStats(inF32, mask, maskVal, stats)) {
     1611            psError(__func__, "p_psVectorClippedStats() failed.\n");
     1612        }
    16461613    }
    16471614    // ************************************************************************

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

    r2204 r2269  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-10-27 00:57:31 $
     12 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-11-03 03:30:30 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    174174);
    175175
    176 
    177 void p_psNormalizeVector(psVector* myData);
    178 
    179 void p_psNormalizeVectorF64(psVector* myData);
     176void psNormalizeVectorRange(psVector* myData,
     177                            float low,
     178                            float high);
    180179
    181180/// @}
Note: See TracChangeset for help on using the changeset viewer.