IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2217


Ignore:
Timestamp:
Oct 27, 2004, 11:06:30 AM (22 years ago)
Author:
gusciora
Message:

Added parameter checking to psFunctions.c

Location:
trunk/psLib/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r2216 r2217  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-27 20:32:49 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-27 21:06:30 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828
    2929/*****************************************************************************
     30 
     31*****************************************************************************/
     32#define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
     33if (NAME < 0) { \
     34    psError(__func__,"%s is less than 0.", #NAME); \
     35    return(RVAL); \
     36}
     37
     38
     39
     40/*****************************************************************************
    3041Macros which take a generic psLib type and determine if it is NULL, or has
    3142the wrong type.
    … …  
    4859    PS_VECTOR macros:
    4960 *****************************************************************************/
    50 #define PS_VECTOR_CHECK_TYPE_F32_OR_F64(VEC, RVAL) \
    51 if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
    52     psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
     61#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
     62if (NAME == NULL || NAME->data.V == NULL) { \
     63    psError(__func__,"Unallowable operation: psVector %s or its data is NULL.", #NAME); \
     64    return(RVAL); \
     65} \
     66
     67#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
     68if (NAME->n < 1) { \
     69    psError(__func__,"Unallowable operation: psVector %s has no elements.", #NAME); \
     70    return(RVAL); \
     71} \
     72
     73#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(NAME, RVAL) \
     74if ((NAME->type.type != PS_TYPE_F32) && (NAME->type.type != PS_TYPE_F64)) { \
     75    psAbort(__func__, "psVector %s: bad type(%d)", #NAME, NAME->type.type); \
    5376    return(RVAL); \
    5477} \
    … …  
    5679#define PS_VECTOR_CHECK_TYPE(NAME, TYPE, RVAL) \
    5780if (NAME->type.type != TYPE) { \
    58     psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
    59     return(RVAL); \
    60 }
    61 
    62 #define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
    63 if (NAME == NULL || NAME->data.V == NULL) { \
    64     psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    65     return(RVAL); \
    66 } \
    67 
    68 #define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
    69 if (NAME->n < 1) { \
    70     psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    71     return(RVAL); \
    72 } \
     81    psError(__func__,"Unallowable operation: psVector %s has incorrect type.", #NAME); \
     82    return(RVAL); \
     83}
    7384
    7485#define PS_VECTOR_CHECK_SIZE_EQUAL(VEC1, VEC2, RVAL) \
    7586if (VEC1->n != VEC2->n) { \
    76     psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     87    psError(__func__,"psVector %s has size %d, psVector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
    7788    return(RVAL); \
    7889}
    … …  
    166177    PS_POLY macros:
    167178*****************************************************************************/
    168 #define PS_POLY_CHECK_NULL(NAME, NULL) \
     179#define PS_POLY_CHECK_NULL(NAME, RVAL) \
    169180if (NAME == NULL || NAME->coeff == NULL) { \
    170181    psError(__func__,"Unallowable operation: polynomial %s or its coeffs is NULL.", #NAME); \
    171     return(NULL); \
     182    return(RVAL); \
    172183} \
    173184

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

    r2204 r2217  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-27 00:57:31 $
     9 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-27 21:06:30 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     13 *
     14 *  XXX: What happens if the polyEVal functions are called with dat of the wrong
     15 *       type?
     16 *
    1317 */
    1418/*****************************************************************************/
    … …  
    3034#include "psAbort.h"
    3135#include "psLogMsg.h"
    32 
    3336#include "psFunctions.h"
     37#include "psConstants.h"
    3438/*****************************************************************************/
    3539/* DEFINE STATEMENTS                                                         */
    3640/*****************************************************************************/
    37 /*
    38  * The following macros take the vector IN as input.  If IN is not of
    39  * the correct type, then a new vector OUT is created of the correct
    40  * type, and the elements of IN and converted to OUT.
    41  *
    42  * If vector IN is the correct type, then OUT is simply set to IN.
    43  */
    44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \
    45 if (IN->type.type != PS_TYPE_F32) { \
    46     psLogMsg(__func__, PS_LOG_WARN, \
    47              "Input vector has incorrect type (%d).\n", \
    48              IN->type); \
    49     OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
    50     \
    51     if (IN->type.type == PS_TYPE_F64) { \
    52         for (psS32 i=0;i<IN->n;i++) { \
    53             OUT->data.F32[i] = (float) (IN->data.F64[i]); \
    54         } \
    55     }\
    56     else { \
    57         psAbort(__func__, "Wrong type.\n"); \
    58     } \
    59 } else { \
    60     OUT = (psVector *) IN; \
    61 }
    62 
    63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \
    64 if (IN->type.type != PS_TYPE_F64) { \
    65     psLogMsg(__func__, PS_LOG_WARN, \
    66              "Input vector has incorrect type (%d).\n", \
    67              IN->type); \
    68     OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
    69     \
    70     if (IN->type.type == PS_TYPE_F32) { \
    71         for (psS32 i=0;i<IN->n;i++) { \
    72             OUT->data.F64[i] = (float) (IN->data.F32[i]); \
    73         } \
    74     }\
    75     else { \
    76         psAbort(__func__, "Wrong type.\n"); \
    77     } \
    78 } else { \
    79     OUT = (psVector *) IN; \
    80 }
    8141
    8242/*****************************************************************************/
    … …  
    11979static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
    12080{
     81    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
     82
    12183    psPolynomial1D **chebPolys = NULL;
    12284    psS32 i = 0;
    … …  
    187149psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
    188150{
     151    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     152
    189153    psVector* gauss = NULL;
    190154    const gsl_rng_type *T = NULL;
    191155    gsl_rng *r = NULL;
    192156    psS32 i = 0;
     157
    193158
    194159    gauss = psVectorAlloc(Npts, PS_TYPE_F32);
    … …  
    212177                                    psPolynomialType type)
    213178{
     179    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     180
    214181    psS32 i = 0;
    215182    psPolynomial1D* newPoly = NULL;
    … …  
    235202                                    psPolynomialType type)
    236203{
     204    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     205    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     206
    237207    psS32 x = 0;
    238208    psS32 y = 0;
    … …  
    268238                                    psPolynomialType type)
    269239{
     240    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     241    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     242    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     243
    270244    psS32 x = 0;
    271245    psS32 y = 0;
    … …  
    310284                                    psPolynomialType type)
    311285{
     286    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     287    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     288    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     289    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     290
    312291    psS32 w = 0;
    313292    psS32 x = 0;
    … …  
    522501float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
    523502{
     503    PS_POLY_CHECK_NULL(myPoly, NAN);
     504
    524505    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    525506        return(p_psOrdPolynomial1DEval(x, myPoly));
    … …  
    535516                                   const psPolynomial1D *myPoly)
    536517{
     518    PS_POLY_CHECK_NULL(myPoly, NULL);
     519    PS_VECTOR_CHECK_NULL(x, NULL);
     520    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     521
    537522    psVector *tmp;
    538523    psVector *myX;
    539524    psS32 i;
    540525
    541     PS_CONVERT_VECTOR_F32(x, myX);
    542 
    543526    tmp = psVectorAlloc(x->n, PS_TYPE_F32);
    544527    for (i=0;i<x->n;i++) {
    … …  
    555538float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    556539{
     540    PS_POLY_CHECK_NULL(myPoly, NAN);
     541
    557542    psS32 loop_x = 0;
    558543    psS32 loop_y = 0;
    … …  
    577562float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    578563{
     564    PS_POLY_CHECK_NULL(myPoly, NAN);
     565
    579566    psS32 loop_x = 0;
    580567    psS32 loop_y = 0;
    … …  
    610597float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    611598{
     599    PS_POLY_CHECK_NULL(myPoly, NAN);
     600
    612601    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    613602        return(p_psOrdPolynomial2DEval(x, y, myPoly));
    … …  
    625614                                   const psPolynomial2D *myPoly)
    626615{
     616    PS_POLY_CHECK_NULL(myPoly, NULL);
     617    PS_VECTOR_CHECK_NULL(x, NULL);
     618    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     619    PS_VECTOR_CHECK_NULL(y, NULL);
     620    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     621
    627622    psVector *tmp;
    628623    psVector *myX;
    … …  
    631626    psS32 vecLen=x->n;
    632627
    633     PS_CONVERT_VECTOR_F32(x, myX);
    634     PS_CONVERT_VECTOR_F32(y, myY);
    635628    vecLen=x->n;
    636629    if (y->n < vecLen) {
    … …  
    727720float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    728721{
     722    PS_POLY_CHECK_NULL(myPoly, NAN);
     723
    729724    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    730725        return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    742737                                   const psPolynomial3D *myPoly)
    743738{
     739    PS_POLY_CHECK_NULL(myPoly, NULL);
     740    PS_VECTOR_CHECK_NULL(x, NULL);
     741    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     742    PS_VECTOR_CHECK_NULL(y, NULL);
     743    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     744    PS_VECTOR_CHECK_NULL(z, NULL);
     745    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     746
    744747    psVector *tmp;
    745748    psVector *myX;
    … …  
    749752    psS32 vecLen=x->n;
    750753
    751     PS_CONVERT_VECTOR_F32(x, myX);
    752     PS_CONVERT_VECTOR_F32(y, myY);
    753     PS_CONVERT_VECTOR_F32(z, myZ);
    754754    vecLen=x->n;
    755755    if (y->n < vecLen) {
    … …  
    869869float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    870870{
     871    PS_POLY_CHECK_NULL(myPoly, NAN);
     872
    871873    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    872874        return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    885887                                   const psPolynomial4D *myPoly)
    886888{
     889    PS_POLY_CHECK_NULL(myPoly, NULL);
     890    PS_VECTOR_CHECK_NULL(w, NULL);
     891    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
     892    PS_VECTOR_CHECK_NULL(x, NULL);
     893    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     894    PS_VECTOR_CHECK_NULL(y, NULL);
     895    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     896    PS_VECTOR_CHECK_NULL(z, NULL);
     897    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     898
    887899    psVector *tmp;
    888900    psVector *myW;
    … …  
    893905    psS32 vecLen=x->n;
    894906
    895     PS_CONVERT_VECTOR_F32(w, myW);
    896     PS_CONVERT_VECTOR_F32(x, myX);
    897     PS_CONVERT_VECTOR_F32(y, myY);
    898     PS_CONVERT_VECTOR_F32(z, myZ);
    899907    vecLen=w->n;
    900908    if (y->n < vecLen) {
    … …  
    939947                                      psPolynomialType type)
    940948{
     949    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     950
    941951    psS32 i = 0;
    942952    psDPolynomial1D* newPoly = NULL;
    … …  
    962972                                      psPolynomialType type)
    963973{
     974    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     975    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     976
    964977    psS32 x = 0;
    965978    psS32 y = 0;
    … …  
    9951008                                      psPolynomialType type)
    9961009{
     1010    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1011    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1012    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1013
    9971014    psS32 x = 0;
    9981015    psS32 y = 0;
    … …  
    10371054                                      psPolynomialType type)
    10381055{
     1056    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     1057    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1058    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1059    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1060
    10391061    psS32 w = 0;
    10401062    psS32 x = 0;
    … …  
    12071229double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
    12081230{
     1231    PS_POLY_CHECK_NULL(myPoly, NAN);
     1232
    12091233    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    12101234        return(p_psDOrdPolynomial1DEval(x, myPoly));
    … …  
    12201244                                    const psDPolynomial1D *myPoly)
    12211245{
     1246    PS_POLY_CHECK_NULL(myPoly, NULL);
     1247    PS_VECTOR_CHECK_NULL(x, NULL);
     1248    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1249
    12221250    psVector *tmp;
    12231251    psVector *myX;
    12241252    psS32 i;
    1225 
    1226     PS_CONVERT_VECTOR_F64(x, myX);
    12271253
    12281254    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    … …  
    12971323double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    12981324{
     1325    PS_POLY_CHECK_NULL(myPoly, NAN);
     1326
    12991327    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    13001328        return(p_psDOrdPolynomial2DEval(x, y, myPoly));
    … …  
    13111339                                    const psDPolynomial2D *myPoly)
    13121340{
     1341    PS_POLY_CHECK_NULL(myPoly, NULL);
     1342    PS_VECTOR_CHECK_NULL(x, NULL);
     1343    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1344    PS_VECTOR_CHECK_NULL(y, NULL);
     1345    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1346
    13131347    psVector *tmp;
    13141348    psVector *myX;
    … …  
    13171351    psS32 vecLen=x->n;
    13181352
    1319     PS_CONVERT_VECTOR_F64(x, myX);
    1320     PS_CONVERT_VECTOR_F64(y, myY);
    13211353    vecLen=x->n;
    13221354    if (y->n < vecLen) {
    … …  
    14131445double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    14141446{
     1447    PS_POLY_CHECK_NULL(myPoly, NAN);
     1448
    14151449    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    14161450        return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    14281462                                    const psDPolynomial3D *myPoly)
    14291463{
     1464    PS_POLY_CHECK_NULL(myPoly, NULL);
     1465    PS_VECTOR_CHECK_NULL(x, NULL);
     1466    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1467    PS_VECTOR_CHECK_NULL(y, NULL);
     1468    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1469    PS_VECTOR_CHECK_NULL(z, NULL);
     1470    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1471
    14301472    psVector *tmp;
    14311473    psVector *myX;
    … …  
    14351477    psS32 vecLen=x->n;
    14361478
    1437     PS_CONVERT_VECTOR_F64(x, myX);
    1438     PS_CONVERT_VECTOR_F64(y, myY);
    1439     PS_CONVERT_VECTOR_F64(z, myZ);
    14401479    vecLen=x->n;
    14411480    if (y->n < vecLen) {
    … …  
    15571596double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    15581597{
     1598    PS_POLY_CHECK_NULL(myPoly, NAN);
     1599
    15591600    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    15601601        return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    15731614                                    const psDPolynomial4D *myPoly)
    15741615{
     1616    PS_POLY_CHECK_NULL(myPoly, NULL);
     1617    PS_VECTOR_CHECK_NULL(w, NULL);
     1618    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
     1619    PS_VECTOR_CHECK_NULL(x, NULL);
     1620    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1621    PS_VECTOR_CHECK_NULL(y, NULL);
     1622    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1623    PS_VECTOR_CHECK_NULL(z, NULL);
     1624    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1625
    15751626    psVector *tmp;
    15761627    psVector *myW;
    … …  
    15811632    psS32 vecLen=x->n;
    15821633
    1583     PS_CONVERT_VECTOR_F64(w, myW);
    1584     PS_CONVERT_VECTOR_F64(x, myX);
    1585     PS_CONVERT_VECTOR_F64(y, myY);
    1586     PS_CONVERT_VECTOR_F64(z, myZ);
    15871634    vecLen=w->n;
    15881635    if (y->n < vecLen) {

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

    r2213 r2217  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-27 20:11:47 $
    13  *
     11 *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-27 21:06:30 $
     13n *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 */
    … …  
    3939/* DEFINE STATEMENTS                                                         */
    4040/*****************************************************************************/
    41 // will use robust statistical methods.
    42 #define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
    43 #define CLIPPED_NUM_ITER_LB 1
    44 #define CLIPPED_NUM_ITER_UB 10
    45 #define CLIPPED_SIGMA_LB 1.0
    46 #define CLIPPED_SIGMA_UB 10.0
    47 #define true 1
    48 #define false 0
    49 #define MY_MAX_FLOAT HUGE
    50 #define MAX_ITERATIONS 10
    51 
    52 void p_psVectorRobustStats(const psVector* restrict myVector,
    53                            const psVector* restrict maskVector, psU32 maskVal, psStats* stats);
    54 
    55 
     41#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
     42#define PS_CLIPPED_NUM_ITER_LB 1
     43#define PS_CLIPPED_NUM_ITER_UB 10
     44#define PS_CLIPPED_SIGMA_LB 1.0
     45#define PS_CLIPPED_SIGMA_UB 10.0
     46#define PS_MAX_FLOAT HUGE
     47#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
     48
     49#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
     50(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
    5651/*****************************************************************************/
    5752/* TYPE DEFINITIONS                                                          */
    … …  
    5954psVector* p_psConvertToF32(psVector* in);
    6055
     56void p_psVectorRobustStats(const psVector* restrict myVector,
     57                           const psVector* restrict maskVector,
     58                           psU32 maskVal,
     59                           psStats* stats);
    6160/*****************************************************************************/
    6261/* GLOBAL VARIABLES                                                          */
    … …  
    245244{
    246245    psS32 i = 0;                  // Loop index variable
    247     float max = -MY_MAX_FLOAT;  // The calculated maximum
     246    float max = -PS_MAX_FLOAT;  // The calculated maximum
    248247    float rangeMin = 0.0;       // Exclude data below this
    249248    float rangeMax = 0.0;       // Exclude date above this
    … …  
    306305{
    307306    psS32 i = 0;                  // Loop index variable
    308     float min = MY_MAX_FLOAT;   // The calculated maximum
     307    float min = PS_MAX_FLOAT;   // The calculated maximum
    309308    float rangeMin = 0.0;       // Exclude data below this
    310309    float rangeMax = 0.0;       // Exclude date above this
    … …  
    493492robustHistogram with a Gaussian of width sigma.
    494493 *****************************************************************************/
    495 #define GAUSS_WIDTH_OTHER 5.0
    496494psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
    497495                                       float sigma)
    … …  
    521519        // warnings?
    522520
    523         x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
     521        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
    524522        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
    525523            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
    … …  
    531529        }
    532530
    533         x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
     531        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
    534532        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
    535533            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
    … …  
    744742
    745743    // Endure that stats->clipIter is within the proper range.
    746     if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
     744    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
    747745        psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
    748746    }
    749747    // Endure that stats->clipSigma is within the proper range.
    750     if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
     748    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
    751749        psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
    752750    }
    … …  
    996994    // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
    997995
    998     while (numIterations < MAX_ITERATIONS) {
     996    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
    999997        midpoint = (rangeHigh + rangeLow) / 2.0;
    1000998        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
    … …  
    10991097    return(tmpFloat);
    11001098}
    1101 
    1102 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
    1103 (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
    11041099
    11051100/******************************************************************************

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

    r2216 r2217  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-27 20:32:49 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-27 21:06:30 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828
    2929/*****************************************************************************
     30 
     31*****************************************************************************/
     32#define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
     33if (NAME < 0) { \
     34    psError(__func__,"%s is less than 0.", #NAME); \
     35    return(RVAL); \
     36}
     37
     38
     39
     40/*****************************************************************************
    3041Macros which take a generic psLib type and determine if it is NULL, or has
    3142the wrong type.
    … …  
    4859    PS_VECTOR macros:
    4960 *****************************************************************************/
    50 #define PS_VECTOR_CHECK_TYPE_F32_OR_F64(VEC, RVAL) \
    51 if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
    52     psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
     61#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
     62if (NAME == NULL || NAME->data.V == NULL) { \
     63    psError(__func__,"Unallowable operation: psVector %s or its data is NULL.", #NAME); \
     64    return(RVAL); \
     65} \
     66
     67#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
     68if (NAME->n < 1) { \
     69    psError(__func__,"Unallowable operation: psVector %s has no elements.", #NAME); \
     70    return(RVAL); \
     71} \
     72
     73#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(NAME, RVAL) \
     74if ((NAME->type.type != PS_TYPE_F32) && (NAME->type.type != PS_TYPE_F64)) { \
     75    psAbort(__func__, "psVector %s: bad type(%d)", #NAME, NAME->type.type); \
    5376    return(RVAL); \
    5477} \
    … …  
    5679#define PS_VECTOR_CHECK_TYPE(NAME, TYPE, RVAL) \
    5780if (NAME->type.type != TYPE) { \
    58     psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
    59     return(RVAL); \
    60 }
    61 
    62 #define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
    63 if (NAME == NULL || NAME->data.V == NULL) { \
    64     psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    65     return(RVAL); \
    66 } \
    67 
    68 #define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
    69 if (NAME->n < 1) { \
    70     psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    71     return(RVAL); \
    72 } \
     81    psError(__func__,"Unallowable operation: psVector %s has incorrect type.", #NAME); \
     82    return(RVAL); \
     83}
    7384
    7485#define PS_VECTOR_CHECK_SIZE_EQUAL(VEC1, VEC2, RVAL) \
    7586if (VEC1->n != VEC2->n) { \
    76     psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     87    psError(__func__,"psVector %s has size %d, psVector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
    7788    return(RVAL); \
    7889}
    … …  
    166177    PS_POLY macros:
    167178*****************************************************************************/
    168 #define PS_POLY_CHECK_NULL(NAME, NULL) \
     179#define PS_POLY_CHECK_NULL(NAME, RVAL) \
    169180if (NAME == NULL || NAME->coeff == NULL) { \
    170181    psError(__func__,"Unallowable operation: polynomial %s or its coeffs is NULL.", #NAME); \
    171     return(NULL); \
     182    return(RVAL); \
    172183} \
    173184

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

    r2204 r2217  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-27 00:57:31 $
     9 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-27 21:06:30 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     13 *
     14 *  XXX: What happens if the polyEVal functions are called with dat of the wrong
     15 *       type?
     16 *
    1317 */
    1418/*****************************************************************************/
    … …  
    3034#include "psAbort.h"
    3135#include "psLogMsg.h"
    32 
    3336#include "psFunctions.h"
     37#include "psConstants.h"
    3438/*****************************************************************************/
    3539/* DEFINE STATEMENTS                                                         */
    3640/*****************************************************************************/
    37 /*
    38  * The following macros take the vector IN as input.  If IN is not of
    39  * the correct type, then a new vector OUT is created of the correct
    40  * type, and the elements of IN and converted to OUT.
    41  *
    42  * If vector IN is the correct type, then OUT is simply set to IN.
    43  */
    44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \
    45 if (IN->type.type != PS_TYPE_F32) { \
    46     psLogMsg(__func__, PS_LOG_WARN, \
    47              "Input vector has incorrect type (%d).\n", \
    48              IN->type); \
    49     OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
    50     \
    51     if (IN->type.type == PS_TYPE_F64) { \
    52         for (psS32 i=0;i<IN->n;i++) { \
    53             OUT->data.F32[i] = (float) (IN->data.F64[i]); \
    54         } \
    55     }\
    56     else { \
    57         psAbort(__func__, "Wrong type.\n"); \
    58     } \
    59 } else { \
    60     OUT = (psVector *) IN; \
    61 }
    62 
    63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \
    64 if (IN->type.type != PS_TYPE_F64) { \
    65     psLogMsg(__func__, PS_LOG_WARN, \
    66              "Input vector has incorrect type (%d).\n", \
    67              IN->type); \
    68     OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
    69     \
    70     if (IN->type.type == PS_TYPE_F32) { \
    71         for (psS32 i=0;i<IN->n;i++) { \
    72             OUT->data.F64[i] = (float) (IN->data.F32[i]); \
    73         } \
    74     }\
    75     else { \
    76         psAbort(__func__, "Wrong type.\n"); \
    77     } \
    78 } else { \
    79     OUT = (psVector *) IN; \
    80 }
    8141
    8242/*****************************************************************************/
    … …  
    11979static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
    12080{
     81    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
     82
    12183    psPolynomial1D **chebPolys = NULL;
    12284    psS32 i = 0;
    … …  
    187149psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
    188150{
     151    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     152
    189153    psVector* gauss = NULL;
    190154    const gsl_rng_type *T = NULL;
    191155    gsl_rng *r = NULL;
    192156    psS32 i = 0;
     157
    193158
    194159    gauss = psVectorAlloc(Npts, PS_TYPE_F32);
    … …  
    212177                                    psPolynomialType type)
    213178{
     179    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     180
    214181    psS32 i = 0;
    215182    psPolynomial1D* newPoly = NULL;
    … …  
    235202                                    psPolynomialType type)
    236203{
     204    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     205    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     206
    237207    psS32 x = 0;
    238208    psS32 y = 0;
    … …  
    268238                                    psPolynomialType type)
    269239{
     240    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     241    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     242    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     243
    270244    psS32 x = 0;
    271245    psS32 y = 0;
    … …  
    310284                                    psPolynomialType type)
    311285{
     286    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     287    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     288    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     289    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     290
    312291    psS32 w = 0;
    313292    psS32 x = 0;
    … …  
    522501float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
    523502{
     503    PS_POLY_CHECK_NULL(myPoly, NAN);
     504
    524505    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    525506        return(p_psOrdPolynomial1DEval(x, myPoly));
    … …  
    535516                                   const psPolynomial1D *myPoly)
    536517{
     518    PS_POLY_CHECK_NULL(myPoly, NULL);
     519    PS_VECTOR_CHECK_NULL(x, NULL);
     520    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     521
    537522    psVector *tmp;
    538523    psVector *myX;
    539524    psS32 i;
    540525
    541     PS_CONVERT_VECTOR_F32(x, myX);
    542 
    543526    tmp = psVectorAlloc(x->n, PS_TYPE_F32);
    544527    for (i=0;i<x->n;i++) {
    … …  
    555538float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    556539{
     540    PS_POLY_CHECK_NULL(myPoly, NAN);
     541
    557542    psS32 loop_x = 0;
    558543    psS32 loop_y = 0;
    … …  
    577562float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    578563{
     564    PS_POLY_CHECK_NULL(myPoly, NAN);
     565
    579566    psS32 loop_x = 0;
    580567    psS32 loop_y = 0;
    … …  
    610597float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    611598{
     599    PS_POLY_CHECK_NULL(myPoly, NAN);
     600
    612601    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    613602        return(p_psOrdPolynomial2DEval(x, y, myPoly));
    … …  
    625614                                   const psPolynomial2D *myPoly)
    626615{
     616    PS_POLY_CHECK_NULL(myPoly, NULL);
     617    PS_VECTOR_CHECK_NULL(x, NULL);
     618    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     619    PS_VECTOR_CHECK_NULL(y, NULL);
     620    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     621
    627622    psVector *tmp;
    628623    psVector *myX;
    … …  
    631626    psS32 vecLen=x->n;
    632627
    633     PS_CONVERT_VECTOR_F32(x, myX);
    634     PS_CONVERT_VECTOR_F32(y, myY);
    635628    vecLen=x->n;
    636629    if (y->n < vecLen) {
    … …  
    727720float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    728721{
     722    PS_POLY_CHECK_NULL(myPoly, NAN);
     723
    729724    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    730725        return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    742737                                   const psPolynomial3D *myPoly)
    743738{
     739    PS_POLY_CHECK_NULL(myPoly, NULL);
     740    PS_VECTOR_CHECK_NULL(x, NULL);
     741    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     742    PS_VECTOR_CHECK_NULL(y, NULL);
     743    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     744    PS_VECTOR_CHECK_NULL(z, NULL);
     745    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     746
    744747    psVector *tmp;
    745748    psVector *myX;
    … …  
    749752    psS32 vecLen=x->n;
    750753
    751     PS_CONVERT_VECTOR_F32(x, myX);
    752     PS_CONVERT_VECTOR_F32(y, myY);
    753     PS_CONVERT_VECTOR_F32(z, myZ);
    754754    vecLen=x->n;
    755755    if (y->n < vecLen) {
    … …  
    869869float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    870870{
     871    PS_POLY_CHECK_NULL(myPoly, NAN);
     872
    871873    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    872874        return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    885887                                   const psPolynomial4D *myPoly)
    886888{
     889    PS_POLY_CHECK_NULL(myPoly, NULL);
     890    PS_VECTOR_CHECK_NULL(w, NULL);
     891    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
     892    PS_VECTOR_CHECK_NULL(x, NULL);
     893    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     894    PS_VECTOR_CHECK_NULL(y, NULL);
     895    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     896    PS_VECTOR_CHECK_NULL(z, NULL);
     897    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     898
    887899    psVector *tmp;
    888900    psVector *myW;
    … …  
    893905    psS32 vecLen=x->n;
    894906
    895     PS_CONVERT_VECTOR_F32(w, myW);
    896     PS_CONVERT_VECTOR_F32(x, myX);
    897     PS_CONVERT_VECTOR_F32(y, myY);
    898     PS_CONVERT_VECTOR_F32(z, myZ);
    899907    vecLen=w->n;
    900908    if (y->n < vecLen) {
    … …  
    939947                                      psPolynomialType type)
    940948{
     949    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     950
    941951    psS32 i = 0;
    942952    psDPolynomial1D* newPoly = NULL;
    … …  
    962972                                      psPolynomialType type)
    963973{
     974    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     975    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     976
    964977    psS32 x = 0;
    965978    psS32 y = 0;
    … …  
    9951008                                      psPolynomialType type)
    9961009{
     1010    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1011    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1012    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1013
    9971014    psS32 x = 0;
    9981015    psS32 y = 0;
    … …  
    10371054                                      psPolynomialType type)
    10381055{
     1056    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     1057    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1058    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1059    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1060
    10391061    psS32 w = 0;
    10401062    psS32 x = 0;
    … …  
    12071229double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
    12081230{
     1231    PS_POLY_CHECK_NULL(myPoly, NAN);
     1232
    12091233    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    12101234        return(p_psDOrdPolynomial1DEval(x, myPoly));
    … …  
    12201244                                    const psDPolynomial1D *myPoly)
    12211245{
     1246    PS_POLY_CHECK_NULL(myPoly, NULL);
     1247    PS_VECTOR_CHECK_NULL(x, NULL);
     1248    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1249
    12221250    psVector *tmp;
    12231251    psVector *myX;
    12241252    psS32 i;
    1225 
    1226     PS_CONVERT_VECTOR_F64(x, myX);
    12271253
    12281254    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    … …  
    12971323double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    12981324{
     1325    PS_POLY_CHECK_NULL(myPoly, NAN);
     1326
    12991327    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    13001328        return(p_psDOrdPolynomial2DEval(x, y, myPoly));
    … …  
    13111339                                    const psDPolynomial2D *myPoly)
    13121340{
     1341    PS_POLY_CHECK_NULL(myPoly, NULL);
     1342    PS_VECTOR_CHECK_NULL(x, NULL);
     1343    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1344    PS_VECTOR_CHECK_NULL(y, NULL);
     1345    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1346
    13131347    psVector *tmp;
    13141348    psVector *myX;
    … …  
    13171351    psS32 vecLen=x->n;
    13181352
    1319     PS_CONVERT_VECTOR_F64(x, myX);
    1320     PS_CONVERT_VECTOR_F64(y, myY);
    13211353    vecLen=x->n;
    13221354    if (y->n < vecLen) {
    … …  
    14131445double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    14141446{
     1447    PS_POLY_CHECK_NULL(myPoly, NAN);
     1448
    14151449    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    14161450        return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    14281462                                    const psDPolynomial3D *myPoly)
    14291463{
     1464    PS_POLY_CHECK_NULL(myPoly, NULL);
     1465    PS_VECTOR_CHECK_NULL(x, NULL);
     1466    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1467    PS_VECTOR_CHECK_NULL(y, NULL);
     1468    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1469    PS_VECTOR_CHECK_NULL(z, NULL);
     1470    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1471
    14301472    psVector *tmp;
    14311473    psVector *myX;
    … …  
    14351477    psS32 vecLen=x->n;
    14361478
    1437     PS_CONVERT_VECTOR_F64(x, myX);
    1438     PS_CONVERT_VECTOR_F64(y, myY);
    1439     PS_CONVERT_VECTOR_F64(z, myZ);
    14401479    vecLen=x->n;
    14411480    if (y->n < vecLen) {
    … …  
    15571596double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    15581597{
     1598    PS_POLY_CHECK_NULL(myPoly, NAN);
     1599
    15591600    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    15601601        return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    15731614                                    const psDPolynomial4D *myPoly)
    15741615{
     1616    PS_POLY_CHECK_NULL(myPoly, NULL);
     1617    PS_VECTOR_CHECK_NULL(w, NULL);
     1618    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
     1619    PS_VECTOR_CHECK_NULL(x, NULL);
     1620    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1621    PS_VECTOR_CHECK_NULL(y, NULL);
     1622    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1623    PS_VECTOR_CHECK_NULL(z, NULL);
     1624    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1625
    15751626    psVector *tmp;
    15761627    psVector *myW;
    … …  
    15811632    psS32 vecLen=x->n;
    15821633
    1583     PS_CONVERT_VECTOR_F64(w, myW);
    1584     PS_CONVERT_VECTOR_F64(x, myX);
    1585     PS_CONVERT_VECTOR_F64(y, myY);
    1586     PS_CONVERT_VECTOR_F64(z, myZ);
    15871634    vecLen=w->n;
    15881635    if (y->n < vecLen) {

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

    r2204 r2217  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-27 00:57:31 $
     9 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-27 21:06:30 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     13 *
     14 *  XXX: What happens if the polyEVal functions are called with dat of the wrong
     15 *       type?
     16 *
    1317 */
    1418/*****************************************************************************/
    … …  
    3034#include "psAbort.h"
    3135#include "psLogMsg.h"
    32 
    3336#include "psFunctions.h"
     37#include "psConstants.h"
    3438/*****************************************************************************/
    3539/* DEFINE STATEMENTS                                                         */
    3640/*****************************************************************************/
    37 /*
    38  * The following macros take the vector IN as input.  If IN is not of
    39  * the correct type, then a new vector OUT is created of the correct
    40  * type, and the elements of IN and converted to OUT.
    41  *
    42  * If vector IN is the correct type, then OUT is simply set to IN.
    43  */
    44 #define PS_CONVERT_VECTOR_F32(IN, OUT) \
    45 if (IN->type.type != PS_TYPE_F32) { \
    46     psLogMsg(__func__, PS_LOG_WARN, \
    47              "Input vector has incorrect type (%d).\n", \
    48              IN->type); \
    49     OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
    50     \
    51     if (IN->type.type == PS_TYPE_F64) { \
    52         for (psS32 i=0;i<IN->n;i++) { \
    53             OUT->data.F32[i] = (float) (IN->data.F64[i]); \
    54         } \
    55     }\
    56     else { \
    57         psAbort(__func__, "Wrong type.\n"); \
    58     } \
    59 } else { \
    60     OUT = (psVector *) IN; \
    61 }
    62 
    63 #define PS_CONVERT_VECTOR_F64(IN, OUT) \
    64 if (IN->type.type != PS_TYPE_F64) { \
    65     psLogMsg(__func__, PS_LOG_WARN, \
    66              "Input vector has incorrect type (%d).\n", \
    67              IN->type); \
    68     OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
    69     \
    70     if (IN->type.type == PS_TYPE_F32) { \
    71         for (psS32 i=0;i<IN->n;i++) { \
    72             OUT->data.F64[i] = (float) (IN->data.F32[i]); \
    73         } \
    74     }\
    75     else { \
    76         psAbort(__func__, "Wrong type.\n"); \
    77     } \
    78 } else { \
    79     OUT = (psVector *) IN; \
    80 }
    8141
    8242/*****************************************************************************/
    … …  
    11979static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
    12080{
     81    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
     82
    12183    psPolynomial1D **chebPolys = NULL;
    12284    psS32 i = 0;
    … …  
    187149psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
    188150{
     151    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     152
    189153    psVector* gauss = NULL;
    190154    const gsl_rng_type *T = NULL;
    191155    gsl_rng *r = NULL;
    192156    psS32 i = 0;
     157
    193158
    194159    gauss = psVectorAlloc(Npts, PS_TYPE_F32);
    … …  
    212177                                    psPolynomialType type)
    213178{
     179    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     180
    214181    psS32 i = 0;
    215182    psPolynomial1D* newPoly = NULL;
    … …  
    235202                                    psPolynomialType type)
    236203{
     204    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     205    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     206
    237207    psS32 x = 0;
    238208    psS32 y = 0;
    … …  
    268238                                    psPolynomialType type)
    269239{
     240    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     241    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     242    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     243
    270244    psS32 x = 0;
    271245    psS32 y = 0;
    … …  
    310284                                    psPolynomialType type)
    311285{
     286    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     287    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     288    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     289    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     290
    312291    psS32 w = 0;
    313292    psS32 x = 0;
    … …  
    522501float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
    523502{
     503    PS_POLY_CHECK_NULL(myPoly, NAN);
     504
    524505    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    525506        return(p_psOrdPolynomial1DEval(x, myPoly));
    … …  
    535516                                   const psPolynomial1D *myPoly)
    536517{
     518    PS_POLY_CHECK_NULL(myPoly, NULL);
     519    PS_VECTOR_CHECK_NULL(x, NULL);
     520    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     521
    537522    psVector *tmp;
    538523    psVector *myX;
    539524    psS32 i;
    540525
    541     PS_CONVERT_VECTOR_F32(x, myX);
    542 
    543526    tmp = psVectorAlloc(x->n, PS_TYPE_F32);
    544527    for (i=0;i<x->n;i++) {
    … …  
    555538float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    556539{
     540    PS_POLY_CHECK_NULL(myPoly, NAN);
     541
    557542    psS32 loop_x = 0;
    558543    psS32 loop_y = 0;
    … …  
    577562float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    578563{
     564    PS_POLY_CHECK_NULL(myPoly, NAN);
     565
    579566    psS32 loop_x = 0;
    580567    psS32 loop_y = 0;
    … …  
    610597float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
    611598{
     599    PS_POLY_CHECK_NULL(myPoly, NAN);
     600
    612601    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    613602        return(p_psOrdPolynomial2DEval(x, y, myPoly));
    … …  
    625614                                   const psPolynomial2D *myPoly)
    626615{
     616    PS_POLY_CHECK_NULL(myPoly, NULL);
     617    PS_VECTOR_CHECK_NULL(x, NULL);
     618    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     619    PS_VECTOR_CHECK_NULL(y, NULL);
     620    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     621
    627622    psVector *tmp;
    628623    psVector *myX;
    … …  
    631626    psS32 vecLen=x->n;
    632627
    633     PS_CONVERT_VECTOR_F32(x, myX);
    634     PS_CONVERT_VECTOR_F32(y, myY);
    635628    vecLen=x->n;
    636629    if (y->n < vecLen) {
    … …  
    727720float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
    728721{
     722    PS_POLY_CHECK_NULL(myPoly, NAN);
     723
    729724    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    730725        return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    742737                                   const psPolynomial3D *myPoly)
    743738{
     739    PS_POLY_CHECK_NULL(myPoly, NULL);
     740    PS_VECTOR_CHECK_NULL(x, NULL);
     741    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     742    PS_VECTOR_CHECK_NULL(y, NULL);
     743    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     744    PS_VECTOR_CHECK_NULL(z, NULL);
     745    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     746
    744747    psVector *tmp;
    745748    psVector *myX;
    … …  
    749752    psS32 vecLen=x->n;
    750753
    751     PS_CONVERT_VECTOR_F32(x, myX);
    752     PS_CONVERT_VECTOR_F32(y, myY);
    753     PS_CONVERT_VECTOR_F32(z, myZ);
    754754    vecLen=x->n;
    755755    if (y->n < vecLen) {
    … …  
    869869float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
    870870{
     871    PS_POLY_CHECK_NULL(myPoly, NAN);
     872
    871873    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    872874        return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    885887                                   const psPolynomial4D *myPoly)
    886888{
     889    PS_POLY_CHECK_NULL(myPoly, NULL);
     890    PS_VECTOR_CHECK_NULL(w, NULL);
     891    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
     892    PS_VECTOR_CHECK_NULL(x, NULL);
     893    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     894    PS_VECTOR_CHECK_NULL(y, NULL);
     895    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     896    PS_VECTOR_CHECK_NULL(z, NULL);
     897    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
     898
    887899    psVector *tmp;
    888900    psVector *myW;
    … …  
    893905    psS32 vecLen=x->n;
    894906
    895     PS_CONVERT_VECTOR_F32(w, myW);
    896     PS_CONVERT_VECTOR_F32(x, myX);
    897     PS_CONVERT_VECTOR_F32(y, myY);
    898     PS_CONVERT_VECTOR_F32(z, myZ);
    899907    vecLen=w->n;
    900908    if (y->n < vecLen) {
    … …  
    939947                                      psPolynomialType type)
    940948{
     949    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
     950
    941951    psS32 i = 0;
    942952    psDPolynomial1D* newPoly = NULL;
    … …  
    962972                                      psPolynomialType type)
    963973{
     974    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     975    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     976
    964977    psS32 x = 0;
    965978    psS32 y = 0;
    … …  
    9951008                                      psPolynomialType type)
    9961009{
     1010    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1011    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1012    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1013
    9971014    psS32 x = 0;
    9981015    psS32 y = 0;
    … …  
    10371054                                      psPolynomialType type)
    10381055{
     1056    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
     1057    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
     1058    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
     1059    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
     1060
    10391061    psS32 w = 0;
    10401062    psS32 x = 0;
    … …  
    12071229double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
    12081230{
     1231    PS_POLY_CHECK_NULL(myPoly, NAN);
     1232
    12091233    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    12101234        return(p_psDOrdPolynomial1DEval(x, myPoly));
    … …  
    12201244                                    const psDPolynomial1D *myPoly)
    12211245{
     1246    PS_POLY_CHECK_NULL(myPoly, NULL);
     1247    PS_VECTOR_CHECK_NULL(x, NULL);
     1248    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1249
    12221250    psVector *tmp;
    12231251    psVector *myX;
    12241252    psS32 i;
    1225 
    1226     PS_CONVERT_VECTOR_F64(x, myX);
    12271253
    12281254    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    … …  
    12971323double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
    12981324{
     1325    PS_POLY_CHECK_NULL(myPoly, NAN);
     1326
    12991327    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    13001328        return(p_psDOrdPolynomial2DEval(x, y, myPoly));
    … …  
    13111339                                    const psDPolynomial2D *myPoly)
    13121340{
     1341    PS_POLY_CHECK_NULL(myPoly, NULL);
     1342    PS_VECTOR_CHECK_NULL(x, NULL);
     1343    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1344    PS_VECTOR_CHECK_NULL(y, NULL);
     1345    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1346
    13131347    psVector *tmp;
    13141348    psVector *myX;
    … …  
    13171351    psS32 vecLen=x->n;
    13181352
    1319     PS_CONVERT_VECTOR_F64(x, myX);
    1320     PS_CONVERT_VECTOR_F64(y, myY);
    13211353    vecLen=x->n;
    13221354    if (y->n < vecLen) {
    … …  
    14131445double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
    14141446{
     1447    PS_POLY_CHECK_NULL(myPoly, NAN);
     1448
    14151449    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    14161450        return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
    … …  
    14281462                                    const psDPolynomial3D *myPoly)
    14291463{
     1464    PS_POLY_CHECK_NULL(myPoly, NULL);
     1465    PS_VECTOR_CHECK_NULL(x, NULL);
     1466    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1467    PS_VECTOR_CHECK_NULL(y, NULL);
     1468    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1469    PS_VECTOR_CHECK_NULL(z, NULL);
     1470    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1471
    14301472    psVector *tmp;
    14311473    psVector *myX;
    … …  
    14351477    psS32 vecLen=x->n;
    14361478
    1437     PS_CONVERT_VECTOR_F64(x, myX);
    1438     PS_CONVERT_VECTOR_F64(y, myY);
    1439     PS_CONVERT_VECTOR_F64(z, myZ);
    14401479    vecLen=x->n;
    14411480    if (y->n < vecLen) {
    … …  
    15571596double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
    15581597{
     1598    PS_POLY_CHECK_NULL(myPoly, NAN);
     1599
    15591600    if (myPoly->type == PS_POLYNOMIAL_ORD) {
    15601601        return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
    … …  
    15731614                                    const psDPolynomial4D *myPoly)
    15741615{
     1616    PS_POLY_CHECK_NULL(myPoly, NULL);
     1617    PS_VECTOR_CHECK_NULL(w, NULL);
     1618    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
     1619    PS_VECTOR_CHECK_NULL(x, NULL);
     1620    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     1621    PS_VECTOR_CHECK_NULL(y, NULL);
     1622    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     1623    PS_VECTOR_CHECK_NULL(z, NULL);
     1624    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
     1625
    15751626    psVector *tmp;
    15761627    psVector *myW;
    … …  
    15811632    psS32 vecLen=x->n;
    15821633
    1583     PS_CONVERT_VECTOR_F64(w, myW);
    1584     PS_CONVERT_VECTOR_F64(x, myX);
    1585     PS_CONVERT_VECTOR_F64(y, myY);
    1586     PS_CONVERT_VECTOR_F64(z, myZ);
    15871634    vecLen=w->n;
    15881635    if (y->n < vecLen) {

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

    r2213 r2217  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-27 20:11:47 $
    13  *
     11 *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-27 21:06:30 $
     13n *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 */
    … …  
    3939/* DEFINE STATEMENTS                                                         */
    4040/*****************************************************************************/
    41 // will use robust statistical methods.
    42 #define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
    43 #define CLIPPED_NUM_ITER_LB 1
    44 #define CLIPPED_NUM_ITER_UB 10
    45 #define CLIPPED_SIGMA_LB 1.0
    46 #define CLIPPED_SIGMA_UB 10.0
    47 #define true 1
    48 #define false 0
    49 #define MY_MAX_FLOAT HUGE
    50 #define MAX_ITERATIONS 10
    51 
    52 void p_psVectorRobustStats(const psVector* restrict myVector,
    53                            const psVector* restrict maskVector, psU32 maskVal, psStats* stats);
    54 
    55 
     41#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
     42#define PS_CLIPPED_NUM_ITER_LB 1
     43#define PS_CLIPPED_NUM_ITER_UB 10
     44#define PS_CLIPPED_SIGMA_LB 1.0
     45#define PS_CLIPPED_SIGMA_UB 10.0
     46#define PS_MAX_FLOAT HUGE
     47#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
     48
     49#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
     50(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
    5651/*****************************************************************************/
    5752/* TYPE DEFINITIONS                                                          */
    … …  
    5954psVector* p_psConvertToF32(psVector* in);
    6055
     56void p_psVectorRobustStats(const psVector* restrict myVector,
     57                           const psVector* restrict maskVector,
     58                           psU32 maskVal,
     59                           psStats* stats);
    6160/*****************************************************************************/
    6261/* GLOBAL VARIABLES                                                          */
    … …  
    245244{
    246245    psS32 i = 0;                  // Loop index variable
    247     float max = -MY_MAX_FLOAT;  // The calculated maximum
     246    float max = -PS_MAX_FLOAT;  // The calculated maximum
    248247    float rangeMin = 0.0;       // Exclude data below this
    249248    float rangeMax = 0.0;       // Exclude date above this
    … …  
    306305{
    307306    psS32 i = 0;                  // Loop index variable
    308     float min = MY_MAX_FLOAT;   // The calculated maximum
     307    float min = PS_MAX_FLOAT;   // The calculated maximum
    309308    float rangeMin = 0.0;       // Exclude data below this
    310309    float rangeMax = 0.0;       // Exclude date above this
    … …  
    493492robustHistogram with a Gaussian of width sigma.
    494493 *****************************************************************************/
    495 #define GAUSS_WIDTH_OTHER 5.0
    496494psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
    497495                                       float sigma)
    … …  
    521519        // warnings?
    522520
    523         x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
     521        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
    524522        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
    525523            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
    … …  
    531529        }
    532530
    533         x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
     531        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
    534532        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
    535533            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
    … …  
    744742
    745743    // Endure that stats->clipIter is within the proper range.
    746     if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
     744    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
    747745        psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
    748746    }
    749747    // Endure that stats->clipSigma is within the proper range.
    750     if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
     748    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
    751749        psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
    752750    }
    … …  
    996994    // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
    997995
    998     while (numIterations < MAX_ITERATIONS) {
     996    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
    999997        midpoint = (rangeHigh + rangeLow) / 2.0;
    1000998        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
    … …  
    10991097    return(tmpFloat);
    11001098}
    1101 
    1102 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
    1103 (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
    11041099
    11051100/******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.