IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6305


Ignore:
Timestamp:
Feb 2, 2006, 11:09:08 AM (20 years ago)
Author:
gusciora
Message:

This file contains some of the generic math functions which
were previously in psStats and psSpline and psMinimize.

Location:
trunk/psLib/src/math
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/Makefile.am

    r6101 r6305  
    1515        psPolynomial.c \
    1616        psSpline.c \
    17         psStats.c
     17        psStats.c \
     18        psMathUtils.c
    1819
    1920EXTRA_DIST = math.i
    … …  
    3233        psPolynomial.h \
    3334        psSpline.h \
    34         psStats.h
     35        psStats.h \
     36        psMathUtils.h

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

    r6226 r6305  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2006-01-27 20:08:58 $
     12 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-02-02 21:09:07 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3434#include "psBinaryOp.h"
    3535#include "psLogMsg.h"
     36#include "psMathUtils.h"
    3637/*****************************************************************************/
    3738/* DEFINE STATEMENTS                                                         */
    … …  
    684685            f->data.F64[i] = y->data.F64[x->n-1];
    685686        } else {
    686             fScalar = p_psVectorInterpolate((psVector *) x, (psVector *) y,
    687                                             3, &tmpScalar);
     687            fScalar = p_psVectorInterpolate(NULL, (psVector *) x, (psVector *) y, 3, &tmpScalar);
    688688            if (fScalar == NULL) {
    689689                psError(PS_ERR_UNKNOWN, false, "Could not perform vector interpolation.  Returning NULL.\n");

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

    r6204 r6305  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-01-26 21:10:22 $
     9*  @version $Revision: 1.141 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-02-02 21:09:07 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    286286
    287287        for (i=nTerms-3;i>=1;i--) {
    288             d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
    289                              (d->data.F64[i+2]);
     288            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
    290289            if(poly->mask[i] == 0) {
    291290                d->data.F64[i] += poly->coeff[i];
    … …  
    293292        }
    294293
    295         tmp = (x * d->data.F64[1]) -
    296               (d->data.F64[2]);
     294        tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
    297295        if(poly->mask[0] == 0) {
    298296            tmp += (0.5 * poly->coeff[0]);

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

    r6204 r6305  
    44*         and evaluation routines.
    55*
    6 *  This file will hold the functions for allocated, freeing, and evaluating
    7 *  splines.
     6*  This file contains the routines that allocate, free, and evaluate splines.
    87*
    9 *  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-01-26 21:10:22 $
     8*  @version $Revision: 1.135 $ $Name: not supported by cvs2svn $
     9*  @date $Date: 2006-02-02 21:09:08 $
    1110*
    1211*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3130#include "psConstants.h"
    3231#include "psErrorText.h"
     32#include "psMathUtils.h"
    3333
    3434/*****************************************************************************/
    … …  
    7272
    7373    return;
     74}
     75
     76
     77static void PS_POLY1D_PRINT(
     78    psPolynomial1D *poly)
     79{
     80    printf("-------------- PS_POLY1D_PRINT() --------------\n");
     81    printf("poly->nX is %d\n", poly->nX);
     82    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
     83        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
     84    }
     85}
     86
     87static void PS_PRINT_SPLINE2(psSpline1D *mySpline)
     88{
     89    printf("-------------- PS_PRINT_SPLINE2() --------------\n");
     90    printf("mySpline->n is %d\n", mySpline->n);
     91    for (psS32 i = 0 ; i < mySpline->n ; i++) {
     92        PS_POLY1D_PRINT(mySpline->spline[i]);
     93    }
     94    PS_VECTOR_PRINT_F32(mySpline->knots);
    7495}
    7596
    … …  
    139160    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    140161    return(derivs2);
    141 }
    142 
    143 
    144 /*****************************************************************************
    145 vectorBinDisectTYPE(): This is a macro for a private function which takes as
    146 input a vector an array of data as well as a single value for that data.  The
    147 input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all
    148 i).  This routine does a binary disection of the vector and returns "i" such
    149 that (v[i] <= x <= v[i+1).  If x lies outside the range of v[], then this
    150 routine prints a warning message and returns (-2 or -1).
    151  *****************************************************************************/
    152 #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
    153 static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
    154                                    psS32 numBins, \
    155                                    ps##TYPE x) \
    156 { \
    157     psS32 min; \
    158     psS32 max; \
    159     psS32 mid; \
    160     psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \
    161     if (x < bins[0]) { \
    162         psLogMsg(__func__, PS_LOG_WARN, \
    163                  "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
    164                  #TYPE, x, bins[0], bins[numBins-1]); \
    165         return(-2); \
    166     } \
    167     if (x > bins[numBins-1]) { \
    168         psLogMsg(__func__, PS_LOG_WARN, \
    169                  "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
    170                  #TYPE, x, bins[0], bins[numBins-1]); \
    171         return(-1); \
    172     } \
    173     \
    174     min = 0; \
    175     max = numBins-2; \
    176     mid = ((max+1)-min)/2; \
    177     while (min != max) { \
    178         psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); \
    179         \
    180         if (x == bins[mid]) { \
    181             psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, mid); \
    182             return(mid); \
    183         } else if (x < bins[mid]) { \
    184             max = mid-1; \
    185         } else { \
    186             min = mid; \
    187         } \
    188         mid = ((max+1)+min)/2; \
    189     } \
    190     psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, min); \
    191     return(min); \
    192 } \
    193 
    194 FUNC_MACRO_VECTOR_BIN_DISECT(S8)
    195 FUNC_MACRO_VECTOR_BIN_DISECT(S16)
    196 FUNC_MACRO_VECTOR_BIN_DISECT(S32)
    197 FUNC_MACRO_VECTOR_BIN_DISECT(S64)
    198 FUNC_MACRO_VECTOR_BIN_DISECT(U8)
    199 FUNC_MACRO_VECTOR_BIN_DISECT(U16)
    200 FUNC_MACRO_VECTOR_BIN_DISECT(U32)
    201 FUNC_MACRO_VECTOR_BIN_DISECT(U64)
    202 FUNC_MACRO_VECTOR_BIN_DISECT(F32)
    203 FUNC_MACRO_VECTOR_BIN_DISECT(F64)
    204 
    205 /*****************************************************************************
    206 p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
    207   *****************************************************************************/
    208 psS32 p_psVectorBinDisect(
    209     psVector *bins,
    210     psScalar *x)
    211 {
    212     PS_ASSERT_VECTOR_NON_NULL(bins, -4);
    213     PS_ASSERT_VECTOR_NON_EMPTY(bins, -4);
    214     PS_ASSERT_PTR_NON_NULL(x, -6);
    215     PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);
    216     char* strType;
    217 
    218     switch (x->type.type) {
    219     case PS_TYPE_U8:
    220         return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));
    221     case PS_TYPE_U16:
    222         return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));
    223     case PS_TYPE_U32:
    224         return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));
    225     case PS_TYPE_U64:
    226         return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));
    227     case PS_TYPE_S8:
    228         return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));
    229     case PS_TYPE_S16:
    230         return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));
    231     case PS_TYPE_S32:
    232         return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
    233     case PS_TYPE_S64:
    234         return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));
    235     case PS_TYPE_F32:
    236         return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
    237     case PS_TYPE_F64:
    238         return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));
    239     case PS_TYPE_C32:
    240         PS_TYPE_NAME(strType,x->type.type);
    241         psError(PS_ERR_BAD_PARAMETER_TYPE,
    242                 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
    243                 strType);
    244         return 0;
    245     case PS_TYPE_C64:
    246         PS_TYPE_NAME(strType,x->type.type);
    247         psError(PS_ERR_BAD_PARAMETER_TYPE,
    248                 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
    249                 strType);
    250         return 0;
    251     case PS_TYPE_BOOL:
    252         PS_TYPE_NAME(strType,x->type.type);
    253         psError(PS_ERR_BAD_PARAMETER_TYPE,
    254                 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
    255                 strType);
    256         return 0;
    257     }
    258     return(-3);
    259 }
    260 
    261 /*****************************************************************************
    262 fullInterpolate1DF32(): This routine will take as input n-element floating
    263 point arrays domain and range, and the x value, assumed to lie with the
    264 domain vector.  It produces as output the (n-1)-order LaGrange interpolated
    265 value of x.
    266  
    267 XXX: do we error check for non-distinct domain values?
    268  
    269 XXX: This is a somewhat general function.  There's no reason to put it in this
    270 file.
    271  
    272 XXX: Sort all these interpolation functions out.  Why are there so many?
    273  
    274 XXX: The fullInterpolate1D macros are not used anywhere.
    275  *****************************************************************************/
    276 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
    277 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
    278                                      ps##TYPE *range, \
    279                                      unsigned int n, \
    280                                      ps##TYPE x) \
    281 { \
    282     \
    283     unsigned int i; \
    284     unsigned int m; \
    285     static psVector *p = NULL; \
    286     p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \
    287     p_psMemSetPersistent(p, true); \
    288     p_psMemSetPersistent(p->data.TYPE, true); \
    289     \
    290     psTrace(__func__, 4, "---- %s() begin %u-order at x=%f) (%d data points) ----\n", __func__, n-1, x, n); \
    291     for (i=0;i<n;i++) { \
    292         psTrace(__func__, 6, "domain/range is (%f %f)\n", domain[i], range[i]); \
    293     } \
    294     \
    295     for (i=0;i<n;i++) { \
    296         p->data.TYPE[i] = range[i]; \
    297         psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \
    298         \
    299     } \
    300     \
    301     /* From NR, during each iteration of the m loop, we are computing the \
    302        p_{i ... i+m} terms. \
    303     */ \
    304     for (m=1;m<n;m++) { \
    305         for (i=0;i<n-m;i++) { \
    306             /* From NR: we are computing P_{i ... i+m} \
    307              */ \
    308             p->data.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \
    309                                ((domain[i]-x) * p->data.TYPE[i+1])) / \
    310                               (domain[i] - domain[i+m]); \
    311             psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \
    312         } \
    313     } \
    314     psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \
    315     return(p->data.TYPE[0]); \
    316 } \
    317 
    318 /* XXX: Do this correctly.
    319 FUNC_MACRO_FULL_INTERPOLATE_1D(U8)
    320 FUNC_MACRO_FULL_INTERPOLATE_1D(U16)
    321 FUNC_MACRO_FULL_INTERPOLATE_1D(U32)
    322 FUNC_MACRO_FULL_INTERPOLATE_1D(U64)
    323 FUNC_MACRO_FULL_INTERPOLATE_1D(S8)
    324 FUNC_MACRO_FULL_INTERPOLATE_1D(S16)
    325 FUNC_MACRO_FULL_INTERPOLATE_1D(S32)
    326 FUNC_MACRO_FULL_INTERPOLATE_1D(S64)
    327 FUNC_MACRO_FULL_INTERPOLATE_1D(F64)
    328 */
    329 FUNC_MACRO_FULL_INTERPOLATE_1D(F32)
    330 
    331 
    332 /*****************************************************************************
    333 interpolate1DF32(): this is the base 1-D flat memory routine to perform
    334 LaGrange interpolation.
    335  *****************************************************************************/
    336 static psF32 interpolate1DF32(
    337     psF32 *domain,
    338     psF32 *range,
    339     unsigned int n,
    340     unsigned int order,
    341     psF32 x)
    342 {
    343     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    344     PS_ASSERT_PTR_NON_NULL(domain, NAN)
    345     PS_ASSERT_PTR_NON_NULL(range, NAN)
    346     // XXX: Check valid values for n, order, and x?
    347 
    348     psS32 binNum;
    349     psS32 numIntPoints = order+1;
    350     psS32 origin;
    351 
    352     binNum = vectorBinDisectF32(domain, n, x);
    353 
    354     if (0 == numIntPoints%2) {
    355         origin = binNum - ((numIntPoints/2) - 1);
    356     } else {
    357         origin = binNum - (numIntPoints/2);
    358         if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
    359             // x is closer to binNum+1.
    360             origin = 1 + (binNum - (numIntPoints/2));
    361         }
    362     }
    363     if (origin < 0) {
    364         origin = 0;
    365     }
    366     if ((origin + numIntPoints) > n) {
    367         origin = n - numIntPoints;
    368     }
    369 
    370     psTrace(__func__, 4, "---- %s(....) end ----\n", __func__);
    371     return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x));
    372 }
    373 
    374 
    375 /*****************************************************************************
    376 p_psVectorInterpolate(): This routine will take as input psVectors domain and
    377 range, and the x value, assumed to lie with the domain vector.  It produces
    378 as output the LaGrange interpolated value of a polynomial of the specified
    379 order around the point x.
    380  
    381 XXX: This stuff does not currently work with a mask.
    382  
    383 XXX: add another psScalar argument for the result (so you don't have to
    384      allocate it each time).
    385  
    386 XXX: The VectorCopy routines seg fault when I declare range32 as static.
    387  *****************************************************************************/
    388 psScalar *p_psVectorInterpolate(
    389     psVector *domain,
    390     psVector *range,
    391     psS32 order,
    392     psScalar *x)
    393 {
    394     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    395     PS_ASSERT_VECTOR_NON_NULL(domain, NULL);
    396     PS_ASSERT_VECTOR_NON_NULL(range, NULL);
    397     PS_ASSERT_PTR_NON_NULL(x, NULL);
    398     PS_ASSERT_INT_NONNEGATIVE(order, NULL);
    399     PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL);
    400     PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);
    401     PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);
    402 
    403     psVector *range32 = NULL;
    404     psVector *domain32 = NULL;
    405     if (order > (domain->n - 1)) {
    406         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    407                 PS_ERRORTEXT_psSpline_NOT_ENOUGH_DATAPOINTS,
    408                 order);
    409         return(NULL);
    410     }
    411 
    412     if (x->type.type == PS_TYPE_F32) {
    413         psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
    414         return(psScalarAlloc(interpolate1DF32(domain->data.F32,
    415                                               range->data.F32,
    416                                               domain->n,
    417                                               order,
    418                                               x->data.F32), PS_TYPE_F32));
    419     } else if (x->type.type == PS_TYPE_F64) {
    420         // XXX: use recycled vectors here.
    421         range32 = psVectorCopy(range32, range, PS_TYPE_F32);
    422         domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
    423 
    424         psScalar *tmpScalar = psScalarAlloc((psF64)
    425                                             interpolate1DF32(domain32->data.F32,
    426                                                              range32->data.F32,
    427                                                              domain32->n,
    428                                                              order,
    429                                                              (psF32) x->data.F64), PS_TYPE_F64);
    430         psFree(range32);
    431         psFree(domain32);
    432 
    433         psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    434         // XXX: Convert data type to F64?
    435         return(tmpScalar);
    436 
    437     } else {
    438         char* strType;
    439         PS_TYPE_NAME(strType,x->type.type);
    440         psError(PS_ERR_BAD_PARAMETER_TYPE,
    441                 PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
    442                 strType);
    443     }
    444 
    445     psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
    446     return(NULL);
    447162}
    448163
    … …  
    627342}
    628343
    629 void PS_POLY1D_PRINT(
    630     psPolynomial1D *poly)
    631 {
    632     printf("-------------- PS_POLY1D_PRINT() --------------\n");
    633     printf("poly->nX is %d\n", poly->nX);
    634     for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
    635         printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
    636     }
    637 }
    638 
    639 void PS_PRINT_SPLINE2(psSpline1D *mySpline)
    640 {
    641     printf("-------------- PS_PRINT_SPLINE2() --------------\n");
    642     printf("mySpline->n is %d\n", mySpline->n);
    643     for (psS32 i = 0 ; i < mySpline->n ; i++) {
    644         PS_POLY1D_PRINT(mySpline->spline[i]);
    645     }
    646     PS_VECTOR_PRINT_F32(mySpline->knots);
    647 }
    648 
    649 
    650 
    651 
    652 
    653 
    654 
    655 
    656 
    657 
    658 
    659 
    660 
    661 
    662 
    663 
    664 
    665 
    666 
    667 
    668 
    669 
    670 
    671 
    672 
    673 
    674 
    675 
    676 
    677 
    678 
    679 
    680 
    681 
    682 
    683 
    684344
    685345/*****************************************************************************
    … …  
    706366
    707367    psS32 n = spline->n;
    708     psS32 binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x);
     368    psScalar tmpScalar;
     369    tmpScalar.type.type = PS_TYPE_F32;
     370    tmpScalar.data.F32 = x;
     371    psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar);
     372
    709373    if (binNum < 0) {
    710374        //
    … …  
    746410    if (psTraceGetLevel(__func__) >= 6) {
    747411        PS_VECTOR_PRINT_F32(x);
    748         // XXX: print spline
     412        PS_PRINT_SPLINE2((psSpline1D *) spline);
    749413    }
    750414

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

    r5517 r6305  
    1212 *  @author GLG, MHPCC
    1313 *
    14  *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
    15  *  @date $Date: 2005-11-15 20:10:32 $
     14 *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2006-02-02 21:09:08 $
    1616 *
    1717 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    9292);
    9393
    94 /** Performs a binary disection on a given vector.
    95  *  Searches through an array of data for a specified value.
    96  *
    97  *  @return psS32    corresponding index number of specified value
    98  */
    99 psS32 p_psVectorBinDisect(
    100     psVector *bins,                    ///< Array of non-decreasing values
    101     psScalar *x                        ///< Target value to find
    102 );
    103 
    104 /** Interpolates a series of data points for evaluation at a specific coordinate.  Uses a
    105  *  Lagrange interpolation method.
    106  *
    107  *  @return psScalar*    Lagrange interpolation value at given location
    108  */
    109 psScalar *p_psVectorInterpolate(
    110     psVector *domain,                  ///< Domain (x coords) for interpolation
    111     psVector *range,                   ///< Range (y coords) for interpolation
    112     psS32 order,                       ///< Order of interpolation function
    113     psScalar *x                        ///< Location at which to evaluate
    114 );
    115 
    11694/** \} */ // End of MathGroup Functions
    11795

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

    r6270 r6305  
    33 *  @ingroup Stat
    44 *
    5  *  This file will hold the definition of the histogram and stats data
     5 *  This file will hold the definitions of the histogram and stats data
    66 *  structures.  It also contains prototypes for procedures which operate
    77 *  on those data structures.
    … …  
    1010 *
    1111 *  XXX: The following stats members are never used, or set in this code.
    12  *      stats->robustN50
    1312 *      stats->clippedNvalues
    14  *      stats->binsize
    1513 *
    1614 * XXX: Must do
    1715 * nSubsample points
    1816 * use ->min and ->max (PS_STAT_USE_RANGE)
    19  * use ->binsize (PS_STAT_USE_BINSIZE)
    2017 *
    21  *
    22  *
    23  *
    24  *  @version $Revision: 1.164 $ $Name: not supported by cvs2svn $
    25  *  @date $Date: 2006-02-01 02:07:24 $
     18 *  @version $Revision: 1.165 $ $Name: not supported by cvs2svn $
     19 *  @date $Date: 2006-02-02 21:09:08 $
    2620 *
    2721 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4943#include "psPolynomial.h"
    5044#include "psConstants.h"
     45#include "psMathUtils.h"
    5146
    5247#include "psErrorText.h"
    … …  
    6156#define PS_CLIPPED_SIGMA_LB 1.0
    6257#define PS_CLIPPED_SIGMA_UB 10.0
    63 #define PS_POLY_MEDIAN_MAX_ITERATIONS 10
     58#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
    6459
    6560#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
    … …  
    847842
    848843    // Calculate the quartile points exactly.
    849     // XXX: We should probably do a bin-midpoint if the quartile is not an
    850     // integer.
    851844    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
    852845    stats->sampleLQ = sortedVector->data.F32[nValues / 4];
    … …  
    859852}
    860853
    861 /******************************************************************************
    862 p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the
    863 stdev of the input vector.
    864 Inputs
    865     myVector
    866     maskVector
    867     maskVal
    868     stats
    869 Returns
    870     NULL
    871  
    872 XXX: remove this
    873  *****************************************************************************/
    874 void p_psVectorSampleStdevOLD(const psVector* myVector,
    875                               const psVector* errors,
    876                               const psVector* maskVector,
    877                               psU32 maskVal,
    878                               psStats* stats)
    879 {
    880     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    881     psS32 i = 0;                  // Loop index variable
    882     psS32 countInt = 0;           // # of data points being used
    883     psF32 countFloat = 0.0;     // # of data points being used
    884     psF32 mean = 0.0;           // The mean
    885     psF32 diff = 0.0;           // Used in calculating stdev
    886     psF32 sumSquares = 0.0;     // temporary variable
    887     psF32 sumDiffs = 0.0;       // temporary variable
    888 
    889     // This procedure requires the mean.  If it has not been already
    890     // calculated, then call p_psVectorSampleMean()
    891     if (0 != isnan(stats->sampleMean)) {
    892         p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
    893     }
    894     // If the mean is NAN, then generate a warning and set the stdev to NAN.
    895     if (0 != isnan(stats->sampleMean)) {
    896         stats->sampleStdev = NAN;
    897         psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
    898         psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    899         return;
    900     }
    901 
    902     mean = stats->sampleMean;
    903     if (stats->options & PS_STAT_USE_RANGE) {
    904         if (maskVector != NULL) {
    905             for (i = 0; i < myVector->n; i++) {
    906                 if (!(maskVal & maskVector->data.U8[i]) &&
    907                         (stats->min <= myVector->data.F32[i]) &&
    908                         (myVector->data.F32[i] <= stats->max)) {
    909                     diff = myVector->data.F32[i] - mean;
    910                     sumSquares += (diff * diff);
    911                     sumDiffs += diff;
    912                     countInt++;
    913                 }
    914             }
    915         } else {
    916             for (i = 0; i < myVector->n; i++) {
    917                 if ((stats->min <= myVector->data.F32[i]) &&
    918                         (myVector->data.F32[i] <= stats->max)) {
    919                     diff = myVector->data.F32[i] - mean;
    920                     sumSquares += (diff * diff);
    921                     sumDiffs += diff;
    922                     countInt++;
    923                 }
    924             }
    925             countInt = myVector->n;
    926         }
    927     } else {
    928         if (maskVector != NULL) {
    929             for (i = 0; i < myVector->n; i++) {
    930                 if (!(maskVal & maskVector->data.U8[i])) {
    931                     diff = myVector->data.F32[i] - mean;
    932                     sumSquares += (diff * diff);
    933                     sumDiffs += diff;
    934                     countInt++;
    935                 }
    936             }
    937         } else {
    938             for (i = 0; i < myVector->n; i++) {
    939                 diff = myVector->data.F32[i] - mean;
    940                 sumSquares += (diff * diff);
    941                 sumDiffs += diff;
    942                 countInt++;
    943             }
    944             countInt = myVector->n;
    945         }
    946     }
    947     if (countInt == 0) {
    948         stats->sampleStdev = NAN;
    949         psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
    950     } else if (countInt == 1) {
    951         stats->sampleStdev = 0.0;
    952         psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
    953     } else {
    954         countFloat = (psF32)countInt;
    955         stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    956     }
    957     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    958 }
    959854/******************************************************************************
    960855p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
    … …  
    1094989    -2: warning
    1095990 
    1096 XXX: Do we really need to calculate median and mean?
    1097  
    1098991XXX: Use static vectors for tmpMask.
    1099  
    1100 XXX: This has not been tested.
    1101992 *****************************************************************************/
    1102993psS32 p_psVectorClippedStats(const psVector* myVector,
    … …  
    12231114        // Otherwise, use the new results and continue.
    12241115        if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
    1225             // Exit loop.  XXX: throw an error/warning here?
     1116            // Exit loop.
    12261117            iter = stats->clipIter;
    1227             rc = -1;
     1118            psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n");
     1119            psFree(tmpMask);
     1120            return(-1);
    12281121        } else {
    12291122            clippedMean = statsTmp->sampleMean;
    … …  
    12461139}
    12471140
    1248 /*****************************************************************************
    1249 These macros and functions define the following functions:
    1250  
    1251     p_psNormalizeVectorRange(myData, low, high)
    1252  
    1253 That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
    1254 can be of any type.  Arguments low/high will be converted to the appropriate
    1255 type and one of the type-specific functions below will be called:
    1256  
    1257     p_psNormalizeVectorRangeU8(myData, low, high)
    1258     p_psNormalizeVectorRangeU16(myData, low, high)
    1259     p_psNormalizeVectorRangeU32(myData, low, high)
    1260     p_psNormalizeVectorRangeU64(myData, low, high)
    1261     p_psNormalizeVectorRangeS8(myData, low, high)
    1262     p_psNormalizeVectorRangeS16(myData, low, high)
    1263     p_psNormalizeVectorRangeS32(myData, low, high)
    1264     p_psNormalizeVectorRangeS64(myData, low, high)
    1265     p_psNormalizeVectorRangeF32(myData, low, high)
    1266     p_psNormalizeVectorRangeF64(myData, low, high)
    1267  *****************************************************************************/
    1268 #define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
    1269 void p_psNormalizeVectorRange##TYPE(psVector* myData, \
    1270                                     ps##TYPE outLow, \
    1271                                     ps##TYPE outHigh) \
    1272 { \
    1273     ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
    1274     ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
    1275     psS32 i = 0; \
    1276     \
    1277     for (i = 0; i < myData->n; i++) { \
    1278         if (myData->data.TYPE[i] < min) { \
    1279             min = myData->data.TYPE[i]; \
    1280         } \
    1281         if (myData->data.TYPE[i] > max) { \
    1282             max = myData->data.TYPE[i]; \
    1283         } \
    1284     } \
    1285     \
    1286     /* Ensure that max!=min before we divide by (max-min) */ \
    1287     if (max != min) { \
    1288         for (i = 0; i < myData->n; i++) { \
    1289             myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
    1290                                     (outHigh - outLow) / (max - min)); \
    1291         } \
    1292     } else { \
    1293         psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
    1294         for (i = 0; i < myData->n; i++) \
    1295         { \
    1296             \
    1297             myData->data.TYPE[i] = outLow; \
    1298             \
    1299         } \
    1300     } \
    1301 } \
    1302 
    1303 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
    1304 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
    1305 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
    1306 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
    1307 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
    1308 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
    1309 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
    1310 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
    1311 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
    1312 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
    1313 
    1314 void p_psNormalizeVectorRange(psVector* myData,
    1315                               psF64 outLow,
    1316                               psF64 outHigh)
    1317 {
    1318     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    1319     switch (myData->type.type) {
    1320     case PS_TYPE_U8:
    1321         p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
    1322         break;
    1323     case PS_TYPE_U16:
    1324         p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
    1325         break;
    1326     case PS_TYPE_U32:
    1327         p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
    1328         break;
    1329     case PS_TYPE_U64:
    1330         p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
    1331         break;
    1332     case PS_TYPE_S8:
    1333         p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
    1334         break;
    1335     case PS_TYPE_S16:
    1336         p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
    1337         break;
    1338     case PS_TYPE_S32:
    1339         p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
    1340         break;
    1341     case PS_TYPE_S64:
    1342         p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
    1343         break;
    1344     case PS_TYPE_F32:
    1345         p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
    1346         break;
    1347     case PS_TYPE_F64:
    1348         p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
    1349         break;
    1350     case PS_TYPE_C32:
    1351     case PS_TYPE_C64:
    1352     default:
    1353         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1354                 "Unallowable operation: %s has incorrect type.",
    1355                 myData);
    1356         break;
    1357     }
    1358     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    1359 }
     1141
     1142
    13601143
    13611144/******************************************************************************
    … …  
    13711154XXX: Create a 2nd-order polynomial version and solve for X analytically.
    13721155 *****************************************************************************/
    1373 psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
    1374                        psF32 rangeLow,
    1375                        psF32 rangeHigh,
    1376                        psF32 getThisValue)
     1156psF32 p_ps1DPolyMedian(
     1157    psPolynomial1D* myPoly,
     1158    psF32 rangeLow,
     1159    psF32 rangeHigh,
     1160    psF32 getThisValue)
    13771161{
    13781162    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    … …  
    14051189
    14061190        f = psPolynomial1DEval(myPoly, midpoint);
     1191        printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f);
    14071192        if (fabs(f - getThisValue) <= FLT_EPSILON) {
    14081193            return (midpoint);
    … …  
    14191204    return (midpoint);
    14201205}
     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)
     1245static psF32 QuadraticInverse(
     1246    psF32 a,
     1247    psF32 b,
     1248    psF32 c,
     1249    psF32 y)
     1250{
     1251    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
     1252
     1253    psF64 x1 = -b/(2.0*a) + tmp;
     1254    psF64 x2 = -b/(2.0*a) - tmp;
     1255
     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);
     1265}
     1266
    14211267
    14221268/******************************************************************************
    … …  
    14381284arbitrary vectors.  We should probably test that condition.
    14391285*****************************************************************************/
    1440 psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
    1441                                         psVector *yVec,
    1442                                         psS32 binNum,
    1443                                         psF32 yVal)
     1286psF32 fitQuadraticSearchForYThenReturnX(
     1287    psVector *xVec,
     1288    psVector *yVec,
     1289    psS32 binNum,
     1290    psF32 yVal)
    14441291{
    14451292    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    … …  
    14611308    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    14621309    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
    1463     psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    14641310
    14651311    psF32 tmpFloat = 0.0f;
    … …  
    14891335            if (!(y->data.F64[1] <= y->data.F64[2])) {
    14901336                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
    1491                 psFree(myPoly);
    14921337                psFree(x);
    14931338                psFree(y);
    … …  
    15051350            if (!(y->data.F64[1] >= y->data.F64[2])) {
    15061351                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
    1507                 psFree(myPoly);
    15081352                psFree(x);
    15091353                psFree(y);
    … …  
    15251369
    15261370        // Determine the coefficients of the polynomial.
    1527         //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
     1371        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    15281372        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
    15291373        if (myPoly == NULL) {
    … …  
    15311375                    false,
    15321376                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
    1533             psFree(myPoly);
    15341377            psFree(x);
    15351378            psFree(y);
    … …  
    15411384        psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
    15421385        psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
    1543         psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
     1386        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n",
     1387                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
    15441388                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
    15451389                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
    … …  
    15481392        // polynomial, looking for the value x such that f(x) = yVal
    15491393        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]);
    15501395        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);
     1399        psFree(myPoly);
     1400
    15511401
    15521402        if (isnan(tmpFloat)) {
    … …  
    15541404                    false,
    15551405                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
    1556             psFree(myPoly);
    15571406            psFree(x);
    15581407            psFree(y);
    … …  
    15811430
    15821431    psTrace(__func__, 6, "FIT: return %f\n", tmpFloat);
    1583     psFree(myPoly);
    15841432    psFree(x);
    15851433    psFree(y);
    … …  
    15911439
    15921440/******************************************************************************
    1593 XXX: We assume unnormalized gaussians.
     1441NOTE: We assume unnormalized gaussians.
    15941442 *****************************************************************************/
    1595 psF32 psMinimizeLMChi2Gauss1D(psVector *deriv,
    1596                               const psVector *params,
    1597                               const psVector *coords)
     1443psF32 psMinimizeLMChi2Gauss1D(
     1444    psVector *deriv,
     1445    const psVector *params,
     1446    const psVector *coords)
    15981447{
    15991448    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    … …  
    16271476}
    16281477
     1478//XXX: Use the psLib function?
    16291479psF32 LinInterpolate(psF32 x0,
    16301480                     psF32 x1,
    … …  
    16431493    return(x0 + y * (x1 - x0) / (y1 - y0));
    16441494}
    1645 
    1646 
    1647 
    1648 
    1649 
    1650 
    1651 
    1652 
    1653 
    1654 
    1655 
    1656 
    1657 
    1658 
    1659 
    1660 
    1661 
    1662 
    1663 
    1664 
    1665 
    1666 
    1667 
    1668 
    1669 
    1670 
    1671 
    1672 
    1673 
    1674 
    1675 
    1676 
    1677 
    1678 
    1679 
    1680 
    1681 
    1682 
    1683 
    1684 
    1685 
    1686 
    1687 
    1688 
    1689 
    1690 
    1691 
    1692 
    1693 
    1694 
    1695 
    1696 
    1697 
    1698 
    1699 
    1700 
    1701 
    1702 
    1703 
    1704 
    1705 
    1706 
    1707 
    1708 
    1709 
    1710 
    1711 
    1712 
    1713 
    1714 
    1715 
    1716 
    17171495
    17181496
    … …  
    17541532XXX: Review and ensure that all memory is free'ed at premature exits.
    17551533*****************************************************************************/
    1756 #define INITIAL_NUM_BINS 1000.0
     1534#define INITIAL_NUM_BINS 500.0
    17571535psS32 p_psVectorRobustStats(const psVector* myVector,
    17581536                            const psVector* errors,
    … …  
    17701548    psHistogram *cumulativeRobustHistogram = NULL;
    17711549    psS32 numBins = 0;
    1772     psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
    1773     tmpScalar->type.type = PS_TYPE_F32;
     1550    psScalar tmpScalar;
     1551    tmpScalar.type.type = PS_TYPE_F32;
    17741552    psS32 totalDataPoints = 0;
    17751553    psS32 rc = 0;
    … …  
    18061584                psFree(tmpStatsMinMax);
    18071585                psFree(tmpMaskVec);
    1808                 psFree(tmpScalar);
    18091586                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    18101587                return(1);
    … …  
    18331610            psFree(tmpStatsMinMax);
    18341611            psFree(tmpMaskVec);
    1835             psFree(tmpScalar);
    18361612
    18371613            psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
    … …  
    18851661            binMedian = 0;
    18861662        } else {
    1887             tmpScalar->data.F32 = totalDataPoints/2.0;
    1888             binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     1663            tmpScalar.data.F32 = totalDataPoints/2.0;
     1664            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
    18891665            if (binMedian < 0) {
    18901666                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
    … …  
    18921668                psFree(robustHistogram);
    18931669                psFree(cumulativeRobustHistogram);
    1894                 psFree(tmpScalar);
    18951670                psFree(tmpMaskVec);
    18961671                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    … …  
    19051680        // XXX: Check return codes.
    19061681        //
     1682        printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n");
    19071683        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
    19081684                                  *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    19181694        psS32 binLo = 0;
    19191695        psS32 binHi = 0;
    1920         tmpScalar->data.F32 = totalDataPoints * 0.158655f;
    1921         if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
     1696        tmpScalar.data.F32 = totalDataPoints * 0.158655f;
     1697        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
    19221698            binLo = 0;
    19231699        } else {
    19241700            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
    1925             binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     1701            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
    19261702            if (binLo > cumulativeRobustHistogram->nums->n-1) {
    19271703                binLo = cumulativeRobustHistogram->nums->n-1;
    19281704            }
    19291705        }
    1930         tmpScalar->data.F32 = totalDataPoints * 0.841345f;
    1931         if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
     1706        tmpScalar.data.F32 = totalDataPoints * 0.841345f;
     1707        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
    19321708            binHi = 0;
    19331709        } else {
    19341710            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
    1935             binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     1711            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
    19361712            if (binHi > cumulativeRobustHistogram->nums->n-1) {
    19371713                binHi = cumulativeRobustHistogram->nums->n-1;
    … …  
    19471723            psFree(robustHistogram);
    19481724            psFree(cumulativeRobustHistogram);
    1949             psFree(tmpScalar);
    19501725            psFree(tmpMaskVec);
    19511726            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    … …  
    20141789        if (sigma < (2 * binSize)) {
    20151790            psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
    2016             psF32 medianLo;
    2017             psF32 medianHi;
    2018             if ((binMedian - 25) > 0) {
    2019                 medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
    2020             } else {
    2021                 medianLo = robustHistogram->bounds->data.F32[0];
    2022             }
    2023             if ((binMedian + 25) < robustHistogram->bounds->n) {
    2024                 medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
    2025             } else {
    2026                 medianHi = robustHistogram->bounds->data.F32[robustHistogram->bounds->n - 1];
    2027             }
    2028 
     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))];
    20291793            psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
    20301794            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);
    … …  
    20521816    psS32 binHi25 = 0;
    20531817
    2054     tmpScalar->data.F32 = totalDataPoints * 0.25f;
    2055     if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
     1818    tmpScalar.data.F32 = totalDataPoints * 0.25f;
     1819    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
    20561820        // XXX: Special case where its in first bin.  Must code last bin.
    20571821        binLo25 = 0;
    20581822    } else {
    2059         binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
    2060     }
    2061     tmpScalar->data.F32 = totalDataPoints * 0.75f;
    2062     if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
     1823        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
     1824    }
     1825    tmpScalar.data.F32 = totalDataPoints * 0.75f;
     1826    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
    20631827        // XXX: Special case where its in first bin.  Must code last bin.
    20641828        binHi25 = 0;
    20651829    } else {
    2066         binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
     1830        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
    20671831    }
    20681832    if ((binLo25 < 0) || (binHi25 < 0)) {
    … …  
    20711835        psFree(robustHistogram);
    20721836        psFree(cumulativeRobustHistogram);
    2073         psFree(tmpScalar);
    20741837        psFree(tmpMaskVec);
    20751838        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    … …  
    20841847    // XXX: Check for errors.
    20851848    //
     1849    printf("ADD: step 8: Interpolate for the lower quartile positions.\n");
    20861850    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
    20871851                           *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    20891853                           binLo25,
    20901854                           totalDataPoints * 0.25f);
     1855    printf("ADD: step 8: Interpolate for the upper quartile positions.\n");
    20911856    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
    20921857                           *(psVector* *)&cumulativeRobustHistogram->bounds,
    … …  
    21331898            psFree(robustHistogram);
    21341899            psFree(cumulativeRobustHistogram);
    2135             psFree(tmpScalar);
    21361900            psFree(tmpMaskVec);
    21371901            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    … …  
    21701934        psS32 binMin = 0;
    21711935        psS32 binMax = 0;
    2172         tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
    2173         if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
     1936        tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma);
     1937        if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) {
    21741938            binMin = 0;
    21751939        } else {
    2176             binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    2177         }
    2178 
    2179         tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
    2180         if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
     1940            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
     1941        }
     1942
     1943        tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma);
     1944        if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
    21811945            binMax = newHistogram->bounds->n-1;
    21821946        } else {
    2183             binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     1947            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
    21841948        }
    21851949        if ((binMin < 0) || (binMax < 0)) {
    … …  
    22071971        // histogram median.
    22081972        //
    2209         tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
    2210         if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
     1973        tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma);
     1974        if (tmpScalar.data.F32 < tmpStatsMinMax->min) {
    22111975            binMin = 0;
    22121976        } else {
    2213             binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     1977            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
    22141978            // XXX: check for errors here.
    22151979        }
    2216         tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
    2217         if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
     1980        tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma);
     1981        if (tmpScalar.data.F32 > tmpStatsMinMax->max) {
    22181982            binMax = newHistogramSmoothed->n - 1;
    22191983        } else {
    2220             binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     1984            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
    22211985            // XXX: check for errors here.
    22221986        }
    … …  
    22752039            psFree(robustHistogram);
    22762040            psFree(cumulativeRobustHistogram);
    2277             psFree(tmpScalar);
    22782041            psFree(newHistogram);
    22792042            psFree(x);
    … …  
    23142077    psFree(tmpStatsMinMax);
    23152078    psFree(cumulativeRobustHistogram);
    2316     psFree(tmpScalar);
    23172079    psFree(tmpMaskVec);
    23182080    psFree(robustHistogram);
    … …  
    23212083    return(0);
    23222084}
    2323 
    2324 
    2325 
    2326 
    2327 
    2328 
    2329 
    2330 
    23312085
    23322086/*****************************************************************************/
    … …  
    23622116    newStruct->robustUQ = NAN;
    23632117    newStruct->robustLQ = NAN;
    2364     newStruct->robustN50 = -1;            // XXX: This is never used
     2118    newStruct->robustN50 = -1;
    23652119    newStruct->fittedMean = NAN;
    23662120    newStruct->fittedStdev = NAN;
    … …  
    23732127    newStruct->min = NAN;
    23742128    newStruct->max = NAN;
    2375     newStruct->binsize = NAN;          // XXX: This is never used
     2129    newStruct->binsize = NAN;
    23762130    newStruct->nSubsample = 100000;
    23772131    newStruct->options = options;

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

    r6215 r6305  
    1414 *  @author GLG, MHPCC
    1515 *
    16  *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
    17  *  @date $Date: 2006-01-26 23:49:11 $
     16 *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
     17 *  @date $Date: 2006-02-02 21:09:08 $
    1818 *
    1919 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    202202);
    203203
    204 // XXX: Create a single, generic, version of the vector normalize function.
    205 // XXX: Ask IfA for a public psLib function.
    206 
    207 /** Normalize the range of a vector containing U8 values  */
    208 void p_psNormalizeVectorRangeU8(
    209     psVector* myData,                  ///< Vector containing the U8 values
    210     psU8 low,                          ///< Minimum value
    211     psU8 high                          ///< Maximum value
    212 );
    213 
    214 /** Normalize the range of a vector containing U16 values  */
    215 void p_psNormalizeVectorRangeU16(
    216     psVector* myData,                  ///< Vector containing the U16 values
    217     psU16 low,                         ///< Minimum value
    218     psU16 high                         ///< Maximum value
    219 );
    220 
    221 /** Normalize the range of a vector containing U32 values  */
    222 void p_psNormalizeVectorRangeU32(
    223     psVector* myData,                  ///< Vector containing the U32 values
    224     psU32 low,                         ///< Minimum value
    225     psU32 high                         ///< Maximum value
    226 );
    227 
    228 /** Normalize the range of a vector containing U64 values  */
    229 void p_psNormalizeVectorRangeU64(
    230     psVector* myData,                  ///< Vector containing the U64 values
    231     psU64 low,                         ///< Minimum value
    232     psU64 high                         ///< Maximum value
    233 );
    234 
    235 /** Normalize the range of a vector containing S8 values  */
    236 void p_psNormalizeVectorRangeS8(
    237     psVector* myData,                  ///< Vector containing the S8 values
    238     psS8 low,                          ///< Minimum value
    239     psS8 high                          ///< Maximum value
    240 );
    241 
    242 /** Normalize the range of a vector containing S16 values  */
    243 void p_psNormalizeVectorRangeS16(
    244     psVector* myData,                  ///< Vector containing the S16 values
    245     psS16 low,                         ///< Minimum value
    246     psS16 high                         ///< Maximum value
    247 );
    248 
    249 /** Normalize the range of a vector containing S32 values  */
    250 void p_psNormalizeVectorRangeS32(
    251     psVector* myData,                  ///< Vector containing the S32 values
    252     psS32 low,                         ///< Minimum value
    253     psS32 high                         ///< Maximum value
    254 );
    255 
    256 /** Normalize the range of a vector containing S64 values  */
    257 void p_psNormalizeVectorRangeS64(
    258     psVector* myData,                  ///< Vector containing the S64 values
    259     psS64 low,                         ///< Minimum value
    260     psS64 high                         ///< Maximum value
    261 );
    262 
    263 /** Normalize the range of a vector containing F32 values  */
    264 void p_psNormalizeVectorRangeF32(
    265     psVector* myData,                  ///< Vector containing the F32 values
    266     psF32 low,                         ///< Minimum value
    267     psF32 high                         ///< Maximum value
    268 );
    269 
    270 /** Normalize the range of a vector containing F64 values  */
    271 void p_psNormalizeVectorRangeF64(
    272     psVector* myData,                  ///< Vector containing the F64 values
    273     psF64 low,                         ///< Minimum value
    274     psF64 high                         ///< Maximum value
    275 );
    276 
    277204/// @}
    278205
Note: See TracChangeset for help on using the changeset viewer.