IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 17, 2006, 3:20:51 PM (20 years ago)
Author:
Paul Price
Message:

Cleaning up code in psStats.c and psMathUtils.[ch]

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/math/psMathUtils.c

    r6305 r7132  
    11/** @file psMathUtils.c
    22 *
    3  *  This file contains standard math rotines.
     3 *  This file contains standard math routines.
    44 *
    5  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-02-02 21:09:07 $
     5 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-05-18 01:20:51 $
    77 *
    88 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    4848
    4949/*****************************************************************************
    50 vectorBinDisectTYPE(): This is a macro for a private function which takes as
    51 input an array of data as well as a single value for that data.  The input
    52 vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i).
    53 This routine does a binary disection of the vector and returns "i" such that
    54 (v[i] <= x <= v[i+1).  If x lies outside the range of v[], then this routine
    55 prints a warning message and returns (-2 or -1).
     50This is a macro covering the various types for the below function.
    5651 *****************************************************************************/
    57 #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
    58 static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
    59                                    psS32 numBins, \
    60                                    ps##TYPE x) \
    61 { \
    62     psS32 min; \
    63     psS32 max; \
    64     psS32 mid; \
     52#define VECTOR_BIN_DISECT_CASE(TYPE) \
     53case PS_TYPE_##TYPE: { \
     54    ps##TYPE *bounds = bins->data.TYPE; \
     55    ps##TYPE value = x->data.TYPE; \
     56    long min; \
     57    long max; \
     58    long mid; \
    6559    psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \
    6660    /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\
    67     if (x < bins[0]) { \
     61    if (value < bounds[0]) { \
    6862        psLogMsg(__func__, PS_LOG_WARN, \
    6963                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
    70                  #TYPE, x, bins[0], bins[numBins-1]); \
     64                 #TYPE, value, bounds[0], bounds[numBins-1]); \
    7165        return(-2); \
    7266    } \
    73     if (x > bins[numBins-1]) { \
     67    if (value > bounds[numBins-1]) { \
    7468        psLogMsg(__func__, PS_LOG_WARN, \
    7569                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
    76                  #TYPE, x, bins[0], bins[numBins-1]); \
     70                 #TYPE, value, bounds[0], bounds[numBins-1]); \
    7771        return(-1); \
    7872    } \
     
    8276    mid = ((max+1)-min)/2; \
    8377    while (min != max) { \
    84         /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\
    8578        \
    86         if (x == bins[mid]) { \
    87             /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\
     79        if (value == bounds[mid]) { \
    8880            psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \
    8981            return(mid); \
    90         } else if (x < bins[mid]) { \
     82        } else if (value < bounds[mid]) { \
    9183            max = mid-1; \
    9284        } else { \
     
    9587        mid = ((max+1)+min)/2; \
    9688    } \
    97     /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \
    9889    psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \
    9990    return(min); \
    100 } \
    101 
    102 FUNC_MACRO_VECTOR_BIN_DISECT(S8)
    103 FUNC_MACRO_VECTOR_BIN_DISECT(S16)
    104 FUNC_MACRO_VECTOR_BIN_DISECT(S32)
    105 FUNC_MACRO_VECTOR_BIN_DISECT(S64)
    106 FUNC_MACRO_VECTOR_BIN_DISECT(U8)
    107 FUNC_MACRO_VECTOR_BIN_DISECT(U16)
    108 FUNC_MACRO_VECTOR_BIN_DISECT(U32)
    109 FUNC_MACRO_VECTOR_BIN_DISECT(U64)
    110 FUNC_MACRO_VECTOR_BIN_DISECT(F32)
    111 FUNC_MACRO_VECTOR_BIN_DISECT(F64)
     91}
    11292
    11393/*****************************************************************************
    114 p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
     94p_psVectorBinDisect(): This function takes as input an array of data as well as a single value for that data.
     95The input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i).  This routine does a
     96binary disection of the vector and returns "i" such that (v[i] <= x <= v[i+1).  If x lies outside the range of
     97v[], then this routine prints a warning message and returns (-2 or -1).
    11598  *****************************************************************************/
    11699psS32 p_psVectorBinDisect(
    117     psVector *bins,
    118     psScalar *x)
     100    const psVector *bins,
     101    const psScalar *x)
    119102{
    120103    PS_ASSERT_VECTOR_NON_NULL(bins, -4);
     
    122105    PS_ASSERT_PTR_NON_NULL(x, -6);
    123106    PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);
    124     char* strType;
     107    long numBins = bins->n;             // Number of bins
    125108
    126109    switch (x->type.type) {
    127     case PS_TYPE_U8:
    128         return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));
    129     case PS_TYPE_U16:
    130         return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));
    131     case PS_TYPE_U32:
    132         return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));
    133     case PS_TYPE_U64:
    134         return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));
    135     case PS_TYPE_S8:
    136         return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));
    137     case PS_TYPE_S16:
    138         return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));
    139     case PS_TYPE_S32:
    140         return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
    141     case PS_TYPE_S64:
    142         return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));
    143     case PS_TYPE_F32:
    144         return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
    145     case PS_TYPE_F64:
    146         return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));
    147     default:
    148         PS_TYPE_NAME(strType, x->type.type);
    149         psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
    150         return -1;
     110        VECTOR_BIN_DISECT_CASE(S8);
     111        VECTOR_BIN_DISECT_CASE(S16);
     112        VECTOR_BIN_DISECT_CASE(S32);
     113        VECTOR_BIN_DISECT_CASE(S64);
     114        VECTOR_BIN_DISECT_CASE(U8);
     115        VECTOR_BIN_DISECT_CASE(U16);
     116        VECTOR_BIN_DISECT_CASE(U32);
     117        VECTOR_BIN_DISECT_CASE(U64);
     118        VECTOR_BIN_DISECT_CASE(F32);
     119        VECTOR_BIN_DISECT_CASE(F64);
     120    default: {
     121            char* strType;
     122            PS_TYPE_NAME(strType, x->type.type);
     123            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
     124            return -1;
     125        }
    151126    }
    152127    return(-3);
     
    154129
    155130
    156 
    157 /*****************************************************************************
    158 fullInterpolateTYPE(): This routine will take as input n-element psVectors
    159 domain and range, and the x value, assumed to lie with the domain vector.  It
    160 produces as output the (order-1)-order LaGrange interpolated value of x.
    161  
    162 XXX: do we error check for non-distinct domain values?
    163  
    164 This routine will only be called inside this file.  So, we do not have to
    165 do PS_ASSERTS.
    166  *****************************************************************************/
    167 #define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \
    168 static psScalar *vectorInterpolate##TYPE( \
    169         psScalar *out, \
    170         psVector *domain, \
    171         psVector *range, \
    172         psS32 order, \
    173         psScalar *x) \
    174 { \
     131/*************************************************************************************************************
     132Helper macro for p_psVectorInterpolate: handles each of the types
     133*************************************************************************************************************/
     134#define VECTOR_INTERPOLATE_CASE(TYPE) \
     135case PS_TYPE_##TYPE: { \
    175136    psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \
    176137    if (x->data.TYPE < domain->data.TYPE[0]) { \
     
    213174    psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \
    214175    return(out); \
    215 } \
    216 
    217 FUNC_MACRO_FULL_INTERPOLATE(U8)
    218 FUNC_MACRO_FULL_INTERPOLATE(U16)
    219 FUNC_MACRO_FULL_INTERPOLATE(U32)
    220 FUNC_MACRO_FULL_INTERPOLATE(U64)
    221 FUNC_MACRO_FULL_INTERPOLATE(S8)
    222 FUNC_MACRO_FULL_INTERPOLATE(S16)
    223 FUNC_MACRO_FULL_INTERPOLATE(S32)
    224 FUNC_MACRO_FULL_INTERPOLATE(S64)
    225 FUNC_MACRO_FULL_INTERPOLATE(F32)
    226 FUNC_MACRO_FULL_INTERPOLATE(F64)
     176}
    227177
    228178/*****************************************************************************
     
    236186psScalar *p_psVectorInterpolate(
    237187    psScalar *out,
    238     psVector *domain,
    239     psVector *range,
     188    const psVector *domain,
     189    const psVector *range,
    240190    psS32 order,
    241     psScalar *x)
     191    const psScalar *x)
    242192{
    243193    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     
    249199    PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);
    250200    PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);
    251     if (out == NULL) {
     201    if (!out) {
    252202        out = psScalarAlloc(0, x->type.type);
    253203    } else {
    254204        PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL);
    255205    }
    256     char* strType;
    257206
    258207    switch (x->type.type) {
    259     case PS_TYPE_U8:
    260         return(vectorInterpolateU8(out, domain, range, order, x));
    261     case PS_TYPE_U16:
    262         return(vectorInterpolateU16(out, domain, range, order, x));
    263     case PS_TYPE_U32:
    264         return(vectorInterpolateU32(out, domain, range, order, x));
    265     case PS_TYPE_U64:
    266         return(vectorInterpolateU64(out, domain, range, order, x));
    267     case PS_TYPE_S8:
    268         return(vectorInterpolateS8(out, domain, range, order, x));
    269     case PS_TYPE_S16:
    270         return(vectorInterpolateS16(out, domain, range, order, x));
    271     case PS_TYPE_S32:
    272         return(vectorInterpolateS32(out, domain, range, order, x));
    273     case PS_TYPE_S64:
    274         return(vectorInterpolateS64(out, domain, range, order, x));
    275     case PS_TYPE_F32:
    276         return(vectorInterpolateF32(out, domain, range, order, x));
    277     case PS_TYPE_F64:
    278         return(vectorInterpolateF64(out, domain, range, order, x));
    279     default:
    280         PS_TYPE_NAME(strType, x->type.type);
    281         psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
    282         return NULL;
     208        VECTOR_INTERPOLATE_CASE(U8);
     209        VECTOR_INTERPOLATE_CASE(U16);
     210        VECTOR_INTERPOLATE_CASE(U32);
     211        VECTOR_INTERPOLATE_CASE(U64);
     212        VECTOR_INTERPOLATE_CASE(S8);
     213        VECTOR_INTERPOLATE_CASE(S16);
     214        VECTOR_INTERPOLATE_CASE(S32);
     215        VECTOR_INTERPOLATE_CASE(S64);
     216        VECTOR_INTERPOLATE_CASE(F32);
     217        VECTOR_INTERPOLATE_CASE(F64);
     218    default: {
     219            char* strType;
     220            PS_TYPE_NAME(strType, x->type.type);
     221            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
     222            return NULL;
     223        }
    283224    }
    284225
     
    287228
    288229/*****************************************************************************
    289 These macros and functions define the following functions:
    290  
    291     p_psNormalizeVectorRange(myData, low, high)
    292  
    293 That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
    294 can be of any type.  Arguments low/high will be converted to the appropriate
    295 type and one of the type-specific functions below will be called:
    296  
    297     p_psNormalizeVectorRangeU8(myData, low, high)
    298     p_psNormalizeVectorRangeU16(myData, low, high)
    299     p_psNormalizeVectorRangeU32(myData, low, high)
    300     p_psNormalizeVectorRangeU64(myData, low, high)
    301     p_psNormalizeVectorRangeS8(myData, low, high)
    302     p_psNormalizeVectorRangeS16(myData, low, high)
    303     p_psNormalizeVectorRangeS32(myData, low, high)
    304     p_psNormalizeVectorRangeS64(myData, low, high)
    305     p_psNormalizeVectorRangeF32(myData, low, high)
    306     p_psNormalizeVectorRangeF64(myData, low, high)
     230Helper macro for p_psNormalizeVectorRange: handles each of the types
    307231 *****************************************************************************/
    308 #define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
    309 static void p_psNormalizeVectorRange##TYPE( \
    310         psVector* myData, \
    311         ps##TYPE outLow, \
    312         ps##TYPE outHigh) \
    313 { \
    314     ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
    315     ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
    316     psS32 i = 0; \
    317     \
    318     for (i = 0; i < myData->n; i++) { \
     232#define NORMALIZE_VECTOR_RANGE_CASE(TYPE) \
     233case PS_TYPE_##TYPE: { \
     234    ps##TYPE low = outLow; \
     235    ps##TYPE high = outHigh; \
     236    ps##TYPE min = (ps##TYPE)PS_MAX_##TYPE; \
     237    ps##TYPE max = (ps##TYPE)-PS_MAX_##TYPE; \
     238    \
     239    for (long i = 0; i < myData->n; i++) { \
    319240        if (myData->data.TYPE[i] < min) { \
    320241            min = myData->data.TYPE[i]; \
     
    327248    /* Ensure that max!=min before we divide by (max-min) */ \
    328249    if (max != min) { \
    329         for (i = 0; i < myData->n; i++) { \
    330             myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
    331                                     (outHigh - outLow) / (max - min)); \
     250        for (long i = 0; i < myData->n; i++) { \
     251            myData->data.TYPE[i] = (low + (myData->data.TYPE[i] - min) * \
     252                                    (high - low) / (max - min)); \
    332253        } \
    333254    } else { \
    334255        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
    335         for (i = 0; i < myData->n; i++) \
     256        for (long i = 0; i < myData->n; i++) \
    336257        { \
    337258            \
    338             myData->data.TYPE[i] = outLow; \
     259            myData->data.TYPE[i] = low; \
    339260            \
    340261        } \
    341262    } \
     263    break; \
    342264} \
    343265
    344 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
    345 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
    346 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
    347 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
    348 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
    349 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
    350 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
    351 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
    352 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
    353 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
    354 
    355 psS32 p_psNormalizeVectorRange(psVector* myData,
    356                                psF64 outLow,
    357                                psF64 outHigh)
     266/*************************************************************************************************************
     267p_psNormalizeVectorRange(myData, low, high): this function normalises the vector (myData) to the range
     268low:
     269high.
     270*************************************************************************************************************/
     271bool p_psNormalizeVectorRange(psVector* myData,
     272                              psF64 outLow,
     273                              psF64 outHigh)
    358274{
    359     PS_ASSERT_VECTOR_NON_NULL(myData, -1);
    360     char* strType;
    361 
     275    PS_ASSERT_VECTOR_NON_NULL(myData, false);
    362276    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     277
    363278    switch (myData->type.type) {
    364     case PS_TYPE_U8:
    365         p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
    366         break;
    367     case PS_TYPE_U16:
    368         p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
    369         break;
    370     case PS_TYPE_U32:
    371         p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
    372         break;
    373     case PS_TYPE_U64:
    374         p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
    375         break;
    376     case PS_TYPE_S8:
    377         p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
    378         break;
    379     case PS_TYPE_S16:
    380         p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
    381         break;
    382     case PS_TYPE_S32:
    383         p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
    384         break;
    385     case PS_TYPE_S64:
    386         p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
    387         break;
    388     case PS_TYPE_F32:
    389         p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
    390         break;
    391     case PS_TYPE_F64:
    392         p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
    393         break;
    394     case PS_TYPE_C32:
    395     case PS_TYPE_C64:
    396     default:
    397         PS_TYPE_NAME(strType, myData->type.type);
    398         psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
    399         break;
     279        NORMALIZE_VECTOR_RANGE_CASE(U8);
     280        NORMALIZE_VECTOR_RANGE_CASE(U16);
     281        NORMALIZE_VECTOR_RANGE_CASE(U32);
     282        NORMALIZE_VECTOR_RANGE_CASE(U64);
     283        NORMALIZE_VECTOR_RANGE_CASE(S8);
     284        NORMALIZE_VECTOR_RANGE_CASE(S16);
     285        NORMALIZE_VECTOR_RANGE_CASE(S32);
     286        NORMALIZE_VECTOR_RANGE_CASE(S64);
     287        NORMALIZE_VECTOR_RANGE_CASE(F32);
     288        NORMALIZE_VECTOR_RANGE_CASE(F64);
     289    default: {
     290            char* strType;
     291            PS_TYPE_NAME(strType, myData->type.type);
     292            psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
     293            break;
     294        }
    400295    }
    401296    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    402     return(0);
    403 }
    404 
    405 
     297    return true;
     298}
     299
     300
Note: See TracChangeset for help on using the changeset viewer.