IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2788


Ignore:
Timestamp:
Dec 21, 2004, 7:09:32 PM (22 years ago)
Author:
gusciora
Message:

Changed types to PS types.

Location:
trunk/psLib/src
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/dataManip/psConstants.h

    r2778 r2788  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-12-21 20:42:07 $
     8 *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-12-22 05:09:32 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     12 *
     13 *  XXX: Add parenthesis around all arguments so that these macros can be
     14 *       called with complex expressions.
     15 *
     16 *
     17 *
    1218 */
    1319
  • trunk/psLib/src/dataManip/psFunctions.c

    r2778 r2788  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-12-21 20:42:07 $
     9 *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-12-22 05:09:32 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5555static void dPolynomial4DFree(psDPolynomial4D* myPoly);
    5656static void spline1DFree(psSpline1D *tmpSpline);
    57 static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x);
     57static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
    5858static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
    5959
     
    277277 
    278278 *****************************************************************************/
    279 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     279static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    280280{
    281281    psS32 loop_x = 0;
    282     float polySum = 0.0;
    283     float xSum = 1.0;
     282    psF32 polySum = 0.0;
     283    psF32 xSum = 1.0;
    284284
    285285    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
     
    307307// XXX: How does the mask vector effect Crenshaw's formula?
    308308// XXX: We assume that x is scaled between -1.0 and 1.0;
    309 static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     309static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    310310{
    311311    //    PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0);
     
    314314    psS32 n;
    315315    psS32 i;
    316     float tmp;
     316    psF32 tmp;
    317317
    318318    n = myPoly->n;
     
    337337    psS32 n;
    338338    psS32 i;
    339     float tmp;
     339    psF32 tmp;
    340340    psPolynomial1D **chebPolys = NULL;
    341341
     
    355355}
    356356
    357 static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     357static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    358358{
    359359    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    361361    psS32 loop_x = 0;
    362362    psS32 loop_y = 0;
    363     float polySum = 0.0;
    364     float xSum = 1.0;
    365     float ySum = 1.0;
     363    psF32 polySum = 0.0;
     364    psF32 xSum = 1.0;
     365    psF32 ySum = 1.0;
    366366
    367367    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    379379}
    380380
    381 static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     381static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    382382{
    383383    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    386386    psS32 loop_y = 0;
    387387    psS32 i = 0;
    388     float polySum = 0.0;
     388    psF32 polySum = 0.0;
    389389    psPolynomial1D* *chebPolys = NULL;
    390390    psS32 maxChebyPoly = 0;
     
    414414}
    415415
    416 static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     416static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    417417{
    418418    psS32 loop_x = 0;
    419419    psS32 loop_y = 0;
    420420    psS32 loop_z = 0;
    421     float polySum = 0.0;
    422     float xSum = 1.0;
    423     float ySum = 1.0;
    424     float zSum = 1.0;
     421    psF32 polySum = 0.0;
     422    psF32 xSum = 1.0;
     423    psF32 ySum = 1.0;
     424    psF32 zSum = 1.0;
    425425
    426426    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    442442}
    443443
    444 static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     444static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    445445{
    446446    psS32 loop_x = 0;
     
    448448    psS32 loop_z = 0;
    449449    psS32 i = 0;
    450     float polySum = 0.0;
     450    psF32 polySum = 0.0;
    451451    psPolynomial1D* *chebPolys = NULL;
    452452    psS32 maxChebyPoly = 0;
     
    483483}
    484484
    485 static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     485static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    486486{
    487487    psS32 loop_w = 0;
     
    489489    psS32 loop_y = 0;
    490490    psS32 loop_z = 0;
    491     float polySum = 0.0;
    492     float wSum = 1.0;
    493     float xSum = 1.0;
    494     float ySum = 1.0;
    495     float zSum = 1.0;
     491    psF32 polySum = 0.0;
     492    psF32 wSum = 1.0;
     493    psF32 xSum = 1.0;
     494    psF32 ySum = 1.0;
     495    psF32 zSum = 1.0;
    496496
    497497    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    517517}
    518518
    519 static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     519static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    520520{
    521521    psS32 loop_w = 0;
     
    524524    psS32 loop_z = 0;
    525525    psS32 i = 0;
    526     float polySum = 0.0;
     526    psF32 polySum = 0.0;
    527527    psPolynomial1D* *chebPolys = NULL;
    528528    psS32 maxChebyPoly = 0;
     
    568568    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    569569 *****************************************************************************/
    570 static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     570static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    571571{
    572572    psS32 loop_x = 0;
    573     double polySum = 0.0;
    574     double xSum = 1.0;
     573    psF64 polySum = 0.0;
     574    psF64 xSum = 1.0;
    575575
    576576    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     
    586586// XXX: You can do this without having to psAlloc() vector d.
    587587// XXX: How does the mask vector effect Crenshaw's formula?
    588 static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     588static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    589589{
    590590    psVector *d;
    591591    psS32 n;
    592592    psS32 i;
    593     double tmp;
     593    psF64 tmp;
    594594
    595595    n = myPoly->n;
     
    611611}
    612612
    613 static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     613static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    614614{
    615615    psS32 loop_x = 0;
    616616    psS32 loop_y = 0;
    617     double polySum = 0.0;
    618     double xSum = 1.0;
    619     double ySum = 1.0;
     617    psF64 polySum = 0.0;
     618    psF64 xSum = 1.0;
     619    psF64 ySum = 1.0;
    620620
    621621    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    633633}
    634634
    635 static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     635static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    636636{
    637637    psS32 loop_x = 0;
    638638    psS32 loop_y = 0;
    639639    psS32 i = 0;
    640     double polySum = 0.0;
     640    psF64 polySum = 0.0;
    641641    psPolynomial1D* *chebPolys = NULL;
    642642    psS32 maxChebyPoly = 0;
     
    667667}
    668668
    669 static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     669static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    670670{
    671671    psS32 loop_x = 0;
    672672    psS32 loop_y = 0;
    673673    psS32 loop_z = 0;
    674     double polySum = 0.0;
    675     double xSum = 1.0;
    676     double ySum = 1.0;
    677     double zSum = 1.0;
     674    psF64 polySum = 0.0;
     675    psF64 xSum = 1.0;
     676    psF64 ySum = 1.0;
     677    psF64 zSum = 1.0;
    678678
    679679    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    695695}
    696696
    697 static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     697static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    698698{
    699699    psS32 loop_x = 0;
     
    701701    psS32 loop_z = 0;
    702702    psS32 i = 0;
    703     double polySum = 0.0;
     703    psF64 polySum = 0.0;
    704704    psPolynomial1D* *chebPolys = NULL;
    705705    psS32 maxChebyPoly = 0;
     
    736736}
    737737
    738 static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     738static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    739739{
    740740    psS32 loop_w = 0;
     
    742742    psS32 loop_y = 0;
    743743    psS32 loop_z = 0;
    744     double polySum = 0.0;
    745     double wSum = 1.0;
    746     double xSum = 1.0;
    747     double ySum = 1.0;
    748     double zSum = 1.0;
     744    psF64 polySum = 0.0;
     745    psF64 wSum = 1.0;
     746    psF64 xSum = 1.0;
     747    psF64 ySum = 1.0;
     748    psF64 zSum = 1.0;
    749749
    750750    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    770770}
    771771
    772 static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     772static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    773773{
    774774    psS32 loop_w = 0;
     
    777777    psS32 loop_z = 0;
    778778    psS32 i = 0;
    779     double polySum = 0.0;
     779    psF64 polySum = 0.0;
    780780    psPolynomial1D* *chebPolys = NULL;
    781781    psS32 maxChebyPoly = 0;
     
    828828 *****************************************************************************/
    829829#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
    830 static float fullInterpolate1D##TYPE(ps##TYPE *domain, \
     830static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
    831831                                     ps##TYPE *range, \
    832832                                     psS32 n, \
     
    896896LaGrange interpolation.
    897897 *****************************************************************************/
    898 static float interpolate1DF32(float *domain,
    899                               float *range,
     898static psF32 interpolate1DF32(psF32 *domain,
     899                              psF32 *range,
    900900                              psS32 n,
    901901                              psS32 order,
    902                               float x)
     902                              psF32 x)
    903903{
    904904    psS32 binNum;
     
    941941    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    942942 *****************************************************************************/
    943 float psGaussian(float x, float mean, float sigma, psBool normal)
    944 {
    945     float tmp = 1.0;
     943psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal)
     944{
     945    psF32 tmp = 1.0;
    946946
    947947    psTrace(".psLib.dataManip.psFunctions.psGaussian", 4,
     
    967967XXX: There is no way to seed the random generator.
    968968 *****************************************************************************/
    969 psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts)
     969psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
    970970{
    971971    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     
    10071007    newPoly->type = type;
    10081008    newPoly->n = n;
    1009     newPoly->coeff = (float *)psAlloc(n * sizeof(float));
    1010     newPoly->coeffErr = (float *)psAlloc(n * sizeof(float));
    1011     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1009    newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32));
     1010    newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32));
     1011    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    10121012    for (i = 0; i < n; i++) {
    10131013        newPoly->coeff[i] = 0.0;
     
    10361036    newPoly->nY = nY;
    10371037
    1038     newPoly->coeff = (float **)psAlloc(nX * sizeof(float *));
    1039     newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *));
    1040     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1038    newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1039    newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1040    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    10411041    for (x = 0; x < nX; x++) {
    1042         newPoly->coeff[x] = (float *)psAlloc(nY * sizeof(float));
    1043         newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float));
    1044         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1042        newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1043        newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1044        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    10451045    }
    10461046    for (x = 0; x < nX; x++) {
     
    10751075    newPoly->nZ = nZ;
    10761076
    1077     newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **));
    1078     newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **));
    1079     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1077    newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1078    newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1079    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    10801080    for (x = 0; x < nX; x++) {
    1081         newPoly->coeff[x] = (float **)psAlloc(nY * sizeof(float *));
    1082         newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *));
    1083         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1081        newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1082        newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1083        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    10841084        for (y = 0; y < nY; y++) {
    1085             newPoly->coeff[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1086             newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1087             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1085            newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1086            newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1087            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    10881088        }
    10891089    }
     
    11241124    newPoly->nZ = nZ;
    11251125
    1126     newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***));
    1127     newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***));
    1128     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1126    newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1127    newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1128    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    11291129    for (w = 0; w < nW; w++) {
    1130         newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **));
    1131         newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **));
    1132         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1130        newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1131        newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1132        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    11331133        for (x = 0; x < nX; x++) {
    1134             newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1135             newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1136             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1134            newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1135            newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1136            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    11371137            for (y = 0; y < nY; y++) {
    1138                 newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1139                 newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1140                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1138                newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1139                newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1140                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    11411141            }
    11421142        }
     
    11571157}
    11581158
    1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, float x)
     1159psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x)
    11601160{
    11611161    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    11951195}
    11961196
    1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y)
     1197psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y)
    11981198{
    11991199    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    12501250}
    12511251
    1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z)
     1252psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z)
    12531253{
    12541254    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    13161316}
    13171317
    1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z)
     1318psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z)
    13191319{
    13201320    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    14061406    newPoly->type = type;
    14071407    newPoly->n = n;
    1408     newPoly->coeff = (double *)psAlloc(n * sizeof(double));
    1409     newPoly->coeffErr = (double *)psAlloc(n * sizeof(double));
    1410     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1408    newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64));
     1409    newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64));
     1410    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    14111411    for (i = 0; i < n; i++) {
    14121412        newPoly->coeff[i] = 0.0;
     
    14351435    newPoly->nY = nY;
    14361436
    1437     newPoly->coeff = (double **)psAlloc(nX * sizeof(double *));
    1438     newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *));
    1439     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1437    newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1438    newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1439    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    14401440    for (x = 0; x < nX; x++) {
    1441         newPoly->coeff[x] = (double *)psAlloc(nY * sizeof(double));
    1442         newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double));
    1443         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1441        newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1442        newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1443        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    14441444    }
    14451445    for (x = 0; x < nX; x++) {
     
    14741474    newPoly->nZ = nZ;
    14751475
    1476     newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **));
    1477     newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **));
    1478     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1476    newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1477    newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1478    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    14791479    for (x = 0; x < nX; x++) {
    1480         newPoly->coeff[x] = (double **)psAlloc(nY * sizeof(double *));
    1481         newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *));
    1482         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1480        newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1481        newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1482        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    14831483        for (y = 0; y < nY; y++) {
    1484             newPoly->coeff[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1485             newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1486             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1484            newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1485            newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1486            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    14871487        }
    14881488    }
     
    15231523    newPoly->nZ = nZ;
    15241524
    1525     newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***));
    1526     newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***));
    1527     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1525    newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1526    newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1527    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    15281528    for (w = 0; w < nW; w++) {
    1529         newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **));
    1530         newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **));
    1531         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1529        newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1530        newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1531        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    15321532        for (x = 0; x < nX; x++) {
    1533             newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1534             newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1535             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1533            newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1534            newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1535            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    15361536            for (y = 0; y < nY; y++) {
    1537                 newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1538                 newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1539                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1537                newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1538                newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1539                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    15401540            }
    15411541        }
     
    15571557
    15581558
    1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x)
     1559psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x)
    15601560{
    15611561    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    15981598
    15991599
    1600 double psDPolynomial2DEval(const psDPolynomial2D* myPoly,
    1601                            double x,
    1602                            double y)
     1600psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly,
     1601                          psF64 x,
     1602                          psF64 y)
    16031603{
    16041604    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    16551655
    16561656
    1657 double psDPolynomial3DEval(const psDPolynomial3D* myPoly,
    1658                            double x,
    1659                            double y,
    1660                            double z)
     1657psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly,
     1658                          psF64 x,
     1659                          psF64 y,
     1660                          psF64 z)
    16611661{
    16621662    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    17241724}
    17251725
    1726 double psDPolynomial4DEval(const psDPolynomial4D* myPoly,
    1727                            double w,
    1728                            double x,
    1729                            double y,
    1730                            double z)
     1726psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly,
     1727                          psF64 w,
     1728                          psF64 x,
     1729                          psF64 y,
     1730                          psF64 z)
    17311731{
    17321732    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    18101810//    psS32 n;
    18111811//    psPolynomial1D **spline;
    1812 //    float *p_psDeriv2;
    1813 //    float *domains;
     1812//    psF32 *p_psDeriv2;
     1813//    psF32 *domains;
    18141814//} psSpline1D;
    18151815
     
    18221822psSpline1D *psSpline1DAlloc(psS32 numSplines,
    18231823                            psS32 order,
    1824                             float min,
    1825                             float max)
     1824                            psF32 min,
     1825                            psF32 max)
    18261826{
    18271827    PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL);
     
    18311831    psSpline1D *tmp = NULL;
    18321832    psS32 i;
    1833     float tmpDomain;
    1834     float width;
     1833    psF32 tmpDomain;
     1834    psF32 width;
    18351835
    18361836    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
     
    18451845    tmp->p_psDeriv2 = NULL;
    18461846
    1847     tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    1848     width = (max - min) / ((float) numSplines);
     1847    tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32));
     1848    width = (max - min) / ((psF32) numSplines);
    18491849
    18501850    (tmp->domains)[0] = min;
     
    18841884    }
    18851885
    1886     tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float));
     1886    tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));
    18871887
    18881888    for (i=0;i<bounds->n;i++) {
     
    19771977{
    19781978    PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3);
    1979     char* strType;
     1979    psS8* strType;
    19801980
    19811981    switch (x->type.type) {
     
    20742074        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
    20752075
    2076         psScalar *tmpScalar = psScalarAlloc((double)
     2076        psScalar *tmpScalar = psScalarAlloc((psF64)
    20772077                                            interpolate1DF32(domain32->data.F32,
    20782078                                                             range32->data.F32,
    20792079                                                             domain32->n,
    20802080                                                             order,
    2081                                                              (float) x->data.F64), PS_TYPE_F64);
     2081                                                             (psF32) x->data.F64), PS_TYPE_F64);
    20822082        psFree(range32);
    20832083        psFree(domain32);
     
    20892089
    20902090    } else {
    2091         char* strType;
     2091        psS8* strType;
    20922092        PS_TYPE_NAME(strType,x->type.type);
    20932093        psError(PS_ERR_BAD_PARAMETER_TYPE,
     
    21152115     the spline fit functions require F32 and F64.
    21162116 *****************************************************************************/
    2117 float psSpline1DEval(
    2118     float x,
     2117psF32 psSpline1DEval(
     2118    psF32 x,
    21192119    const psSpline1D *spline
    21202120)
     
    21722172        for (i=0;i<x->n;i++) {
    21732173            tmpVector->data.F32[i] = psSpline1DEval(
    2174                                          (float) x->data.F64[i],
     2174                                         (psF32) x->data.F64[i],
    21752175                                         spline
    21762176                                     );
    21772177        }
    21782178    } else {
    2179         char* strType;
     2179        psS8* strType;
    21802180        PS_TYPE_NAME(strType,x->type.type);
    21812181        psError(PS_ERR_BAD_PARAMETER_TYPE,
  • trunk/psLib/src/dataManip/psFunctions.h

    r2600 r2788  
    1212*  @author GLG, MHPCC
    1313*
    14 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-12-02 21:12:51 $
     14*  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-12-22 05:09:32 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3939 *        \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    4040 *
    41  *  @return float      value on the gaussian curve given the input parameters
    42  */
    43 float psGaussian(
    44     float x,                           ///< Value at which to evaluate
    45     float mean,                        ///< Mean for the Gaussian
    46     float stddev,                      ///< Standard deviation for the Gaussian
     41 *  @return psF32      value on the gaussian curve given the input parameters
     42 */
     43psF32 psGaussian(
     44    psF32 x,                           ///< Value at which to evaluate
     45    psF32 mean,                        ///< Mean for the Gaussian
     46    psF32 stddev,                      ///< Standard deviation for the Gaussian
    4747    psBool normal                        ///< Indicates whether result should be normalized
    4848);
     
    5555 */
    5656psVector* p_psGaussianDev(
    57     float mean,                        ///< The mean of the Gaussian
    58     float sigma,                       ///< The sigma of the Gaussian
     57    psF32 mean,                        ///< The mean of the Gaussian
     58    psF32 sigma,                       ///< The sigma of the Gaussian
    5959    psS32 Npts                           ///< The size of the vector
    6060);
     
    7070    psPolynomialType type;             ///< Polynomial type
    7171    psS32 n;                             ///< Number of terms
    72     float *coeff;                      ///< Coefficients
    73     float *coeffErr;                   ///< Error in coefficients
    74     char *mask;                        ///< Coefficient mask
     72    psF32 *coeff;                      ///< Coefficients
     73    psF32 *coeffErr;                   ///< Error in coefficients
     74    psU8 *mask;                        ///< Coefficient mask
    7575}
    7676psPolynomial1D;
     
    8282    psS32 nX;                            ///< Number of terms in x
    8383    psS32 nY;                            ///< Number of terms in y
    84     float **coeff;                     ///< Coefficients
    85     float **coeffErr;                  ///< Error in coefficients
    86     char **mask;                       ///< Coefficients mask
     84    psF32 **coeff;                     ///< Coefficients
     85    psF32 **coeffErr;                  ///< Error in coefficients
     86    psU8 **mask;                       ///< Coefficients mask
    8787}
    8888psPolynomial2D;
     
    9595    psS32 nY;                            ///< Number of terms in y
    9696    psS32 nZ;                            ///< Number of terms in z
    97     float ***coeff;                    ///< Coefficients
    98     float ***coeffErr;                 ///< Error in coefficients
    99     char ***mask;                      ///< Coefficients mask
     97    psF32 ***coeff;                    ///< Coefficients
     98    psF32 ***coeffErr;                 ///< Error in coefficients
     99    psU8 ***mask;                      ///< Coefficients mask
    100100}
    101101psPolynomial3D;
     
    109109    psS32 nY;                            ///< Number of terms in y
    110110    psS32 nZ;                            ///< Number of terms in z
    111     float ****coeff;                   ///< Coefficients
    112     float ****coeffErr;                ///< Error in coefficients
    113     char ****mask;                     ///< Coefficients mask
     111    psF32 ****coeff;                   ///< Coefficients
     112    psF32 ****coeffErr;                ///< Error in coefficients
     113    psU8 ****mask;                     ///< Coefficients mask
    114114}
    115115psPolynomial4D;
     
    160160/** Evaluates a 1-D polynomial at specific coordinates.
    161161 *
    162  *  @return float    result of polynomial at given location
    163  */
    164 float psPolynomial1DEval(
     162 *  @return psF32    result of polynomial at given location
     163 */
     164psF32 psPolynomial1DEval(
    165165    const psPolynomial1D* myPoly,       ///< Coefficients for the polynomial
    166     float x                           ///< location at which to evaluate
     166    psF32 x                           ///< location at which to evaluate
    167167);
    168168
    169169/** Evaluates a 2-D polynomial at specific coordinates.
    170170 *
    171  *  @return float    result of polynomial at given location
    172  */
    173 float psPolynomial2DEval(
     171 *  @return psF32    result of polynomial at given location
     172 */
     173psF32 psPolynomial2DEval(
    174174    const psPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    175     float x,                           ///< x location at which to evaluate
    176     float y                           ///< y location at which to evaluate
     175    psF32 x,                           ///< x location at which to evaluate
     176    psF32 y                           ///< y location at which to evaluate
    177177);
    178178
    179179/** Evaluates a 3-D polynomial at specific coordinates.
    180180 *
    181  *  @return float    result of polynomial at given location
    182  */
    183 float psPolynomial3DEval(
     181 *  @return psF32    result of polynomial at given location
     182 */
     183psF32 psPolynomial3DEval(
    184184    const psPolynomial3D* myPoly,       ///< Coefficients for the polynomial
    185     float x,                           ///< x location at which to evaluate
    186     float y,                           ///< y location at which to evaluate
    187     float z                           ///< z location at which to evaluate
     185    psF32 x,                           ///< x location at which to evaluate
     186    psF32 y,                           ///< y location at which to evaluate
     187    psF32 z                           ///< z location at which to evaluate
    188188);
    189189
    190190/** Evaluates a 4-D polynomial at specific coordinates.
    191191 *
    192  *  @return float    result of polynomial at given location
    193  */
    194 float psPolynomial4DEval(
     192 *  @return psF32    result of polynomial at given location
     193 */
     194psF32 psPolynomial4DEval(
    195195    const psPolynomial4D* myPoly,       ///< Coefficients for the polynomial
    196     float w,                           ///< w location at which to evaluate
    197     float x,                           ///< x location at which to evaluate
    198     float y,                           ///< y location at which to evaluate
    199     float z                           ///< z location at which to evaluate
     196    psF32 w,                           ///< w location at which to evaluate
     197    psF32 x,                           ///< x location at which to evaluate
     198    psF32 y,                           ///< y location at which to evaluate
     199    psF32 z                           ///< z location at which to evaluate
    200200);
    201201
     
    235235    psPolynomialType type;             ///< Polynomial type
    236236    psS32 n;                             ///< Number of terms
    237     double *coeff;                     ///< Coefficients
    238     double *coeffErr;                  ///< Error in coefficients
    239     char *mask;                        ///< Coefficient mask
     237    psF64 *coeff;                     ///< Coefficients
     238    psF64 *coeffErr;                  ///< Error in coefficients
     239    psU8 *mask;                        ///< Coefficient mask
    240240}
    241241psDPolynomial1D;
     
    247247    psS32 nX;                            ///< Number of terms in x
    248248    psS32 nY;                            ///< Number of terms in y
    249     double **coeff;                    ///< Coefficients
    250     double **coeffErr;                 ///< Error in coefficients
    251     char **mask;                       ///< Coefficients mask
     249    psF64 **coeff;                    ///< Coefficients
     250    psF64 **coeffErr;                 ///< Error in coefficients
     251    psU8 **mask;                       ///< Coefficients mask
    252252}
    253253psDPolynomial2D;
     
    260260    psS32 nY;                            ///< Number of terms in y
    261261    psS32 nZ;                            ///< Number of terms in z
    262     double ***coeff;                   ///< Coefficients
    263     double ***coeffErr;                ///< Error in coefficients
    264     char ***mask;                      ///< Coefficient mask
     262    psF64 ***coeff;                   ///< Coefficients
     263    psF64 ***coeffErr;                ///< Error in coefficients
     264    psU8 ***mask;                      ///< Coefficient mask
    265265}
    266266psDPolynomial3D;
     
    274274    psS32 nY;                            ///< Number of terms in y
    275275    psS32 nZ;                            ///< Number of terms in z
    276     double ****coeff;                  ///< Coefficients
    277     double ****coeffErr;               ///< Error in coefficients
    278     char ****mask;                     ///< Coefficients mask
     276    psF64 ****coeff;                  ///< Coefficients
     277    psF64 ****coeffErr;               ///< Error in coefficients
     278    psU8 ****mask;                     ///< Coefficients mask
    279279}
    280280psDPolynomial4D;
     
    324324/** Evaluates a double-precision 1-D polynomial at specific coordinates.
    325325 *
    326  *  @return float    result of polynomial at given location
    327  */
    328 double psDPolynomial1DEval(
     326 *  @return psF32    result of polynomial at given location
     327 */
     328psF64 psDPolynomial1DEval(
    329329    const psDPolynomial1D* myPoly,      ///< Coefficients for the polynomial
    330     double x                          ///< Value at which to evaluate
     330    psF64 x                          ///< Value at which to evaluate
    331331);
    332332
    333333/** Evaluates a double-precision 2-D polynomial at specific coordinates.
    334334 *
    335  *  @return float    result of polynomial at given location
    336  */
    337 double psDPolynomial2DEval(
     335 *  @return psF32    result of polynomial at given location
     336 */
     337psF64 psDPolynomial2DEval(
    338338    const psDPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    339     double x,                           ///< Value x at which to evaluate
    340     double y            ///< Value y at which to evaluate
     339    psF64 x,                           ///< Value x at which to evaluate
     340    psF64 y            ///< Value y at which to evaluate
    341341);
    342342
    343343/** Evaluates a double-precision 3-D polynomial at specific coordinates.
    344344 *
    345  *  @return float    result of polynomial at given location
    346  */
    347 double psDPolynomial3DEval(
     345 *  @return psF32    result of polynomial at given location
     346 */
     347psF64 psDPolynomial3DEval(
    348348    const psDPolynomial3D* myPoly,      ///< Coefficients for the polynomial
    349     double x,                          ///< Value x at which to evaluate
    350     double y,                          ///< Value y at which to evaluate
    351     double z     ///< Value z at which to evaluate
     349    psF64 x,                          ///< Value x at which to evaluate
     350    psF64 y,                          ///< Value y at which to evaluate
     351    psF64 z     ///< Value z at which to evaluate
    352352);
    353353
    354354/** Evaluates a double-precision 4-D polynomial at specific coordinates.
    355355 *
    356  *  @return float    result of polynomial at given location
    357  */
    358 double psDPolynomial4DEval(
     356 *  @return psF32    result of polynomial at given location
     357 */
     358psF64 psDPolynomial4DEval(
    359359    const psDPolynomial4D* myPoly,      ///< Coefficients for the polynomial
    360     double w,                          ///< Value w at which to evaluate
    361     double x,                          ///< Value x at which to evaluate
    362     double y,                          ///< Value y at which to evaluate
    363     double z     ///< Value z at which to evaluate
     360    psF64 w,                          ///< Value w at which to evaluate
     361    psF64 x,                          ///< Value x at which to evaluate
     362    psF64 y,                          ///< Value y at which to evaluate
     363    psF64 z     ///< Value z at which to evaluate
    364364);
    365365
     
    396396    psS32 n;                        ///< The number of spline polynomials
    397397    psPolynomial1D **spline;      ///< An array of n pointers to the spline polynomials
    398     float *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    399     float *domains;               ///< The boundaries between each spline piece.  Size is n+1.
     398    psF32 *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
     399    psF32 *domains;               ///< The boundaries between each spline piece.  Size is n+1.
    400400}
    401401psSpline1D;
     
    403403psSpline1D *psSpline1DAlloc(psS32 n,
    404404                            psS32 order,
    405                             float min,
    406                             float max);
     405                            psF32 min,
     406                            psF32 max);
    407407
    408408psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
    409409                                   psS32 order);
    410410
    411 float psSpline1DEval(
    412     float x,
     411psF32 psSpline1DEval(
     412    psF32 x,
    413413    const psSpline1D *spline
    414414);
     
    427427                                psScalar *x);
    428428
    429 float p_psNRSpline1DEval(psSpline1D *spline,
     429psF32 p_psNRSpline1DEval(psSpline1D *spline,
    430430                         const psVector* restrict x,
    431431                         const psVector* restrict y,
    432                          float X);
     432                         psF32 X);
    433433
    434434/* \} */// End of MathGroup Functions
  • trunk/psLib/src/dataManip/psMinimize.c

    r2741 r2788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-17 00:18:31 $
     11 *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-22 05:09:32 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 *  XXX: must follow coding name standards on local functions.
    17  *
    18  *  XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via
    19  *       LMM) is not addressed here.  We are waiting for subsequent SDRs
    20  *       which will redefine the LMM functions.
    2117 *
    2218 */
     
    6460XXX: Use a static vector.
    6561 *****************************************************************************/
    66 void psBuildSums1D(double x,
     62void psBuildSums1D(psF64 x,
    6763                   psS32 polyOrder,
    6864                   psVector* sums)
    6965{
    7066    psS32 i = 0;
    71     double xSum = 0.0;
     67    psF64 xSum = 0.0;
    7268
    7369    if (sums == NULL) {
     
    10096XXX: do an F64 version?
    10197 *****************************************************************************/
    102 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
     98psF32 *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    10399                             const psVector* restrict y)        ///< Coordinates
    104100{
     
    108104    psS32 i;
    109105    psS32 k;
    110     float sig;
    111     float p;
     106    psF32 sig;
     107    psF32 p;
    112108    psS32 n = y->n;
    113     float *u = (float *) psAlloc(n * sizeof(float));
    114     float *derivs2 = (float *) psAlloc(n * sizeof(float));
    115     float *X = (float *) & (x->data.F32[0]);
    116     float *Y = (float *) & (y->data.F32[0]);
    117     float qn;
     109    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
     110    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
     111    psF32 *X = (psF32 *) & (x->data.F32[0]);
     112    psF32 *Y = (psF32 *) & (y->data.F32[0]);
     113    psF32 qn;
    118114
    119115    // XXX: The second derivatives at the endpoints, undefined in the SDR,
     
    168164 *****************************************************************************/
    169165/*
    170 float p_psNRSpline1DEval(psSpline1D *spline,
     166psF32 p_psNRSpline1DEval(psSpline1D *spline,
    171167                         const psVector* restrict x,
    172168                         const psVector* restrict y,
    173                          float X)
     169                         psF32 X)
    174170{
    175171    PS_PTR_CHECK_NULL(spline, NAN);
     
    185181    psS32 klo;
    186182    psS32 khi;
    187     float H;
    188     float A;
    189     float B;
    190     float C;
    191     float D;
    192     float Y;
     183    psF32 H;
     184    psF32 A;
     185    psF32 B;
     186    psF32 C;
     187    psF32 D;
     188    psF32 Y;
    193189 
    194190    n = spline->n;
     
    258254            "---- psVectorFitSpline1D() begin ----\n");
    259255    psS32 numSplines = (y->n)-1;
    260     float tmp;
    261     float H;
     256    psF32 tmp;
     257    psF32 H;
    262258    psS32 i;
    263     float slope;
     259    psF32 slope;
    264260    psVector *x32 = NULL;
    265261    psVector *y32 = NULL;
     
    421417    psTrace(".psLib.dataManip.psMinimize", 4,
    422418            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
    423     float x;
    424     int i;
    425     float mean = params->data.F32[0];
    426     float stdev = params->data.F32[1];
     419    psF32 x;
     420    psS32 i;
     421    psF32 mean = params->data.F32[0];
     422    psF32 stdev = params->data.F32[1];
    427423    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    428424
     
    443439    for (i=0;i<coords->n;i++) {
    444440        x = ((psVector *) (coords->data[i]))->data.F32[0];
    445         float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
     441        psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    446442        deriv->data.F32[i][0] = tmp / (stdev * stdev);
    447443        tmp = (x - mean) * (x - mean) *
     
    491487    PS_PTR_CHECK_NULL(params, NULL);
    492488
    493     double normalization = params->data.F32[0];
    494     double x0 = params->data.F32[1];
    495     double y0 = params->data.F32[2];
    496     double sigmaX = params->data.F32[3];
    497     double sigmaY = params->data.F32[4];
    498     double theta = params->data.F32[5];
     489    psF64 normalization = params->data.F32[0];
     490    psF64 x0 = params->data.F32[1];
     491    psF64 y0 = params->data.F32[2];
     492    psF64 sigmaX = params->data.F32[3];
     493    psF64 sigmaY = params->data.F32[4];
     494    psF64 theta = params->data.F32[5];
    499495    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    500496
     
    502498        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
    503499    } else {
    504         // XXX: Check size of derivative
     500        PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL);
    505501    }
    506502
     
    508504            "---- psMinimizeLMChi2Gauss2D() begin ----\n");
    509505
    510     for (int i=0;i<coords->n;i++) {
    511         double x = ((psVector *) coords->data[i])->data.F32[0];
    512         double y = ((psVector *) coords->data[i])->data.F32[0];
    513 
    514         double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
    515         double v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
    516 
    517         double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
    518                                              v*v/(2.0 * sigmaY * sigmaY)))/
    519                       (2.0 * M_PI * sigmaX * sigmaY);
     506    for (psS32 i=0;i<coords->n;i++) {
     507        psF64 x = ((psVector *) coords->data[i])->data.F32[0];
     508        psF64 y = ((psVector *) coords->data[i])->data.F32[0];
     509
     510        psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
     511        psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
     512
     513        psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
     514                                            v*v/(2.0 * sigmaY * sigmaY)))/
     515                     (2.0 * M_PI * sigmaX * sigmaY);
    520516        out->data.F32[i] = flux;
    521517
     
    584580        PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL);
    585581    }
    586 
    587     // XXX: Generate code that modifies covar matrix if not-NULL.
    588582    if (covar != NULL) {
    589583        PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL);
     
    610604    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    611605    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    612 
    613     //    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
    614     //    for (i=0;i<numData;i++) {
    615     //        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
    616     //    }
    617606    psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32);
    618 
    619     psVector *currValueVec;
    620     psVector *newValueVec;
    621 
    622     float currChi2 = 0.0;
    623     float newChi2 = 0.0;
    624     float lamda = 0.00005;
     607    psVector *currValueVec = NULL;
     608    psVector *newValueVec = NULL;
     609    psF32 currChi2 = 0.0;
     610    psF32 newChi2 = 0.0;
     611    psF32 lamda = 0.00005;
    625612    lamda = 0.05;
    626613
     
    711698                if (j == k) {
    712699                    A->data.F64[j][k] =
    713                         (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
     700                        (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]);
    714701                } else {
    715                     A->data.F64[j][k] = (double) alpha->data.F32[j][k];
     702                    A->data.F64[j][k] = (psF64) alpha->data.F32[j][k];
    716703                }
    717704            }
     
    741728            } else {
    742729                newParams->data.F32[i] = params->data.F32[i] -
    743                                          (float) paramDeltasF64->data.F64[i];
     730                                         (psF32) paramDeltasF64->data.F64[i];
    744731            }
    745732        }
     
    777764            // We already masked params.
    778765            for (i=0;i<numParams;i++) {
    779                 params->data.F32[i] = (float) newParams->data.F32[i];
     766                params->data.F32[i] = (psF32) newParams->data.F32[i];
    780767            }
    781768            lamda*= 0.1;
     
    826813    psS32 k;
    827814    psS32 n = x->n;
    828     double fac;
    829     double sum;
     815    psF64 fac;
     816    psF64 sum;
    830817    PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64);
    831818    psScalar *fScalar;
     
    836823    // variable declarations.  I retain them here to maintain coherence
    837824    // with the NR code.
    838     double min = -1.0;
    839     double max = 1.0;
    840     double bma = 0.5 * (max-min);  // 1
    841     double bpa = 0.5 * (max+min);  // 0
     825    psF64 min = -1.0;
     826    psF64 max = 1.0;
     827    psF64 bma = 0.5 * (max-min);  // 1
     828    psF64 bpa = 0.5 * (max+min);  // 0
    842829
    843830    // In this loop, we first calculate the values of X for which the
     
    851838    for (psS32 i=0;i<n;i++) {
    852839        // NR 5.8.4
    853         double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
    854         double X = (Y + bma + bpa) - 1.0;
     840        psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n));
     841        psF64 X = (Y + bma + bpa) - 1.0;
    855842        tmpScalar.data.F64 = X;
    856843
     
    873860    // coefficients of the Chebyshev polynomial: NR 5.8.7.
    874861
    875     fac = 2.0/((float) n);
     862    fac = 2.0/((psF32) n);
    876863    // XXX: is this loop bound correct?
    877864    for (j=0;j<myPoly->n;j++) {
     
    879866        for (k=0;k<n;k++) {
    880867            sum+= f->data.F64[k] *
    881                   cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
     868                  cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n));
    882869        }
    883870
     
    10751062 *****************************************************************************/
    10761063psMinimization *psMinimizationAlloc(psS32 maxIter,
    1077                                     float tol)
     1064                                    psF32 tol)
    10781065{
    10791066    PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL);
     
    10921079// LINE to it.  We assume BASEMASK is non-null.
    10931080#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
    1094 for (int i=0;i<BASE->n;i++) { \
     1081for (psS32 i=0;i<BASE->n;i++) { \
    10951082    if (BASEMASK->data.U8[i] == 0) { \
    10961083        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
     
    11021089#define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \
    11031090BOOL_VAR = true; \
    1104 for (int i=0;i<IN->n;i++) { \
     1091for (psS32 i=0;i<IN->n;i++) { \
    11051092    if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \
    11061093        BOOL_VAR = false; \
     
    11111098#define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \
    11121099BOOL_VAR = true; \
    1113 for (int i=0;i<IN->n;i++) { \
     1100for (psS32 i=0;i<IN->n;i++) { \
    11141101    if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \
    11151102        BOOL_VAR = false; \
     
    11521139                               psMinimizePowellFunc func)
    11531140{
    1154     float a = 0.0;
    1155     float b = 0.0;
    1156     float c = 0.0;
    1157     float fa = 0.0;
    1158     float fb = 0.0;
    1159     float fc = 0.0;
     1141    psF32 a = 0.0;
     1142    psF32 b = 0.0;
     1143    psF32 c = 0.0;
     1144    psF32 fa = 0.0;
     1145    psF32 fb = 0.0;
     1146    psF32 fc = 0.0;
    11601147    psS32 iter = 100;
    1161     float aDir = 0.0;
    1162     float cDir = 0.0;
    1163     float new_aDir = 0.0;
    1164     float new_cDir = 0.0;
     1148    psF32 aDir = 0.0;
     1149    psF32 cDir = 0.0;
     1150    psF32 new_aDir = 0.0;
     1151    psF32 new_cDir = 0.0;
    11651152    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1166     float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1153    psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
    11671154    psVector *tmp = NULL;
    11681155    psBool boolLineIsNull = true;
     
    13011288                                psMinimizePowellFunc func)
    13021289{
    1303     float a = 0.0;
    1304     float b = 0.0;
    1305     float c = 0.0;
    1306     float fa = 0.0;
    1307     float fb = 0.0;
    1308     float fc = 0.0;
     1290    psF32 a = 0.0;
     1291    psF32 b = 0.0;
     1292    psF32 c = 0.0;
     1293    psF32 fa = 0.0;
     1294    psF32 fb = 0.0;
     1295    psF32 fc = 0.0;
    13091296    psS32 iter = 0;
    13101297    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
    13111298    psBool boolLineIsNull = true;
    1312     float prevMin = 0.0;
    1313     int countMin = 0;
     1299    psF32 prevMin = 0.0;
     1300    psS32 countMin = 0;
    13141301
    13151302    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     
    14031390 *****************************************************************************/
    14041391#define PS_LINEMIN_MAX_ITERATIONS 30
    1405 float p_psLineMin(psMinimization *min,
     1392psF32 p_psLineMin(psMinimization *min,
    14061393                  psVector *params,
    14071394                  psVector *line,
     
    14231410    PS_PTR_CHECK_NULL(func, NAN);
    14241411    psVector *bracket;
    1425     float a = 0.0;
    1426     float b = 0.0;
    1427     float c = 0.0;
    1428     float n = 0.0;
    1429     float fa = 0.0;
    1430     float fb = 0.0;
    1431     float fc = 0.0;
    1432     float fn = 0.0;
    1433     float mul = 0.0;
     1412    psF32 a = 0.0;
     1413    psF32 b = 0.0;
     1414    psF32 c = 0.0;
     1415    psF32 n = 0.0;
     1416    psF32 fa = 0.0;
     1417    psF32 fb = 0.0;
     1418    psF32 fc = 0.0;
     1419    psF32 fn = 0.0;
     1420    psF32 mul = 0.0;
    14341421    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
    14351422    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
     
    15441531This routine must minimize a possibly multi-dimensional function.  The
    15451532function to be minimized "func" is:
    1546     float func(psVector *params, psArray *coords)
     1533    psF32 func(psVector *params, psArray *coords)
    15471534The "params" are the parameters of the function which are varied.  The data
    15481535points at which the function is varied are in the argument "coords" which is
     
    15781565    psVector *myParamMask = NULL;
    15791566    psMinimization dummyMin;
    1580     float mul = 0.0;
    1581     float baseFuncVal = 0.0;
    1582     float currFuncVal = 0.0;
     1567    psF32 mul = 0.0;
     1568    psF32 baseFuncVal = 0.0;
     1569    psF32 currFuncVal = 0.0;
    15831570    psS32 biggestIter = 0;
    1584     float biggestDiff = 0.0;
    1585     int iterationNumber = 0;
     1571    psF32 biggestDiff = 0.0;
     1572    psS32 iterationNumber = 0;
    15861573
    15871574    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     
    17111698            }
    17121699        }
    1713         float fqp = func(pQP, coords);
    1714         float term1 = (baseFuncVal - currFuncVal) - biggestDiff;
     1700        psF32 fqp = func(pQP, coords);
     1701        psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff;
    17151702        term1*= term1;
    17161703        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
    1717         float term2 = baseFuncVal - fqp;
     1704        psF32 term2 = baseFuncVal - fqp;
    17181705        term2*= term2 * biggestDiff;
    17191706        if (term1 < term2) {
     
    17601747    PS_PTR_CHECK_NULL(params, NULL);
    17611748
    1762     float x;
    1763     int i;
    1764     float mean = params->data.F32[0];
    1765     float stdev = params->data.F32[1];
     1749    psF32 x;
     1750    psS32 i;
     1751    psF32 mean = params->data.F32[0];
     1752    psF32 stdev = params->data.F32[1];
    17661753    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    17671754
     
    17831770XXX: This is F32 only
    17841771 *****************************************************************************/
    1785 float myPowellChi2Func(const psVector *params,
     1772psF32 myPowellChi2Func(const psVector *params,
    17861773                       const psArray *coords)
    17871774{
     
    17941781    PS_PTR_CHECK_NULL(coords, NAN);
    17951782
    1796     float chi2 = 0.0;
    1797     float d;
     1783    psF32 chi2 = 0.0;
     1784    psF32 d;
    17981785    psS32 i;
    17991786    psVector *tmp;
  • trunk/psLib/src/dataManip/psMinimize.h

    r2250 r2788  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-11-01 23:57:08 $
     10 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-12-22 05:09:32 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4040{
    4141    psS32 maxIter;                       ///< Convergence limit
    42     float tol;                         ///< Error Tolerance
    43     float value;                       ///< Value of function at minimum
     42    psF32 tol;                         ///< Error Tolerance
     43    psF32 value;                       ///< Value of function at minimum
    4444    psS32 iter;                          ///< Number of iterations required
    45     float lastDelta;                   ///< The last difference for the fit
     45    psF32 lastDelta;                   ///< The last difference for the fit
    4646}
    4747psMinimization;
    4848
    4949psMinimization *psMinimizationAlloc(psS32 maxIter,
    50                                     float tol);
     50                                    psF32 tol);
    5151
    5252/** Derive a polynomial fit.
     
    9191
    9292typedef
    93 float (*psMinimizePowellFunc)(const psVector *params,
     93psF32 (*psMinimizePowellFunc)(const psVector *params,
    9494                              const psArray *coords);
    9595
  • trunk/psLib/src/dataManip/psStats.c

    r2782 r2788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.107 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-22 00:54:28 $
     11 *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-22 05:09:32 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    7070/*****************************************************************************/
    7171
    72 psBool p_psGetStatValue(const psStats* stats, double *value)
     72psBool p_psGetStatValue(const psStats* stats, psF64 *value)
    7373{
    7474
     
    145145this routine sets stats->sampleMean to NAN.
    146146 *****************************************************************************/
    147 int p_psVectorSampleMean(const psVector* restrict myVector,
    148                          const psVector* restrict errors,
    149                          const psVector* restrict maskVector,
    150                          psU32 maskVal,
    151                          psStats* stats)
     147psS32 p_psVectorSampleMean(const psVector* restrict myVector,
     148                           const psVector* restrict errors,
     149                           const psVector* restrict maskVector,
     150                           psU32 maskVal,
     151                           psStats* stats)
    152152{
    153153
    154154    psS32 i = 0;                // Loop index variable
    155     float mean = 0.0;           // The mean
     155    psF32 mean = 0.0;           // The mean
    156156    psS32 count = 0;            // # of points in this mean
    157157
     
    171171                }
    172172                if (count != 0) {
    173                     mean /= (float)count;
     173                    mean /= (psF32)count;
    174174                } else {
    175175                    mean = NAN;
     
    185185                }
    186186                if (count != 0) {
    187                     mean /= (float)count;
     187                    mean /= (psF32)count;
    188188                } else {
    189189                    mean = NAN;
     
    199199                }
    200200                if (count != 0) {
    201                     mean /= (float)count;
     201                    mean /= (psF32)count;
    202202                } else {
    203203                    mean = NAN;
     
    207207                    mean += myVector->data.F32[i];
    208208                }
    209                 mean /= (float)myVector->n;
     209                mean /= (psF32)myVector->n;
    210210            }
    211211        }
     
    283283this routine sets stats->max to NAN.
    284284 *****************************************************************************/
    285 int p_psVectorMax(const psVector* restrict myVector,
    286                   const psVector* restrict maskVector,
    287                   psU32 maskVal,
    288                   psStats* stats)
     285psS32 p_psVectorMax(const psVector* restrict myVector,
     286                    const psVector* restrict maskVector,
     287                    psU32 maskVal,
     288                    psStats* stats)
    289289{
    290290    psS32 i = 0;                // Loop index variable
    291     float max = -PS_MAX_F32;    // The calculated maximum
     291    psF32 max = -PS_MAX_F32;    // The calculated maximum
    292292    psS32 empty = true;         // Does this vector have valid elements?
    293293
     
    348348this routine sets stats->min to NAN.
    349349 *****************************************************************************/
    350 int p_psVectorMin(const psVector* restrict myVector,
    351                   const psVector* restrict maskVector,
    352                   psU32 maskVal,
    353                   psStats* stats)
     350psS32 p_psVectorMin(const psVector* restrict myVector,
     351                    const psVector* restrict maskVector,
     352                    psU32 maskVal,
     353                    psStats* stats)
    354354{
    355355    psS32 i = 0;                // Loop index variable
    356     float min = PS_MAX_F32;   // The calculated maximum
     356    psF32 min = PS_MAX_F32;   // The calculated maximum
    357357    psS32 empty = true;         // Does this vector have valid elements?
    358358
     
    612612 *****************************************************************************/
    613613psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
    614                                        float sigma)
     614                                       psF32 sigma)
    615615{
    616616    PS_PTR_CHECK_NULL(robustHistogram, NULL);
     
    619619    psS32 i = 0;                  // Loop index variable
    620620    psS32 j = 0;                  // Loop index variable
    621     float iMid;
    622     float jMid;
     621    psF32 iMid;
     622    psF32 jMid;
    623623    psS32 numBins = robustHistogram->nums->n;
    624624    psS32 numBounds = robustHistogram->bounds->n;
     
    626626    psS32 jMin = 0;
    627627    psS32 jMax = 0;
    628     float firstBound = robustHistogram->bounds->data.F32[0];
    629     float lastBound = robustHistogram->bounds->data.F32[numBounds-1];
     628    psF32 firstBound = robustHistogram->bounds->data.F32[0];
     629    psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1];
    630630    psScalar x;
    631631
     
    795795    psS32 i = 0;                  // Loop index variable
    796796    psS32 countInt = 0;           // # of data points being used
    797     float countFloat = 0.0;     // # of data points being used
    798     float mean = 0.0;           // The mean
    799     float diff = 0.0;           // Used in calculating stdev
    800     float sumSquares = 0.0;     // temporary variable
    801     float sumDiffs = 0.0;       // temporary variable
     797    psF32 countFloat = 0.0;     // # of data points being used
     798    psF32 mean = 0.0;           // The mean
     799    psF32 diff = 0.0;           // Used in calculating stdev
     800    psF32 sumSquares = 0.0;     // temporary variable
     801    psF32 sumDiffs = 0.0;       // temporary variable
    802802
    803803    // This procedure requires the mean.  If it has not been already
     
    865865        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
    866866    } else {
    867         countFloat = (float)countInt;
     867        countFloat = (psF32)countInt;
    868868        stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    869869    }
     
    889889    psS32 i = 0;                  // Loop index variable
    890890    psS32 countInt = 0;           // # of data points being used
    891     float countFloat = 0.0;     // # of data points being used
    892     float mean = 0.0;           // The mean
    893     float diff = 0.0;           // Used in calculating stdev
    894     float sumSquares = 0.0;     // temporary variable
    895     float sumDiffs = 0.0;       // temporary variable
     891    psF32 countFloat = 0.0;     // # of data points being used
     892    psF32 mean = 0.0;           // The mean
     893    psF32 diff = 0.0;           // Used in calculating stdev
     894    psF32 sumSquares = 0.0;     // temporary variable
     895    psF32 sumDiffs = 0.0;       // temporary variable
    896896    //    psF32 sum1;
    897897    //    psF32 sum2;
     
    982982            stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor));
    983983        } else {
    984             countFloat = (float)countInt;
     984            countFloat = (psF32)countInt;
    985985            stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    986986
     
    10031003    -2: warning
    10041004 *****************************************************************************/
    1005 int p_psVectorClippedStats(const psVector* restrict myVector,
    1006                            const psVector* restrict errors,
    1007                            const psVector* restrict maskVector,
    1008                            psU32 maskVal,
    1009                            psStats* stats)
     1005psS32 p_psVectorClippedStats(const psVector* restrict myVector,
     1006                             const psVector* restrict errors,
     1007                             const psVector* restrict maskVector,
     1008                             psU32 maskVal,
     1009                             psStats* stats)
    10101010{
    10111011    psS32 i = 0;                  // Loop index variable
    10121012    psS32 j = 0;                  // Loop index variable
    1013     float clippedMean = 0.0;    // self-explanatory
    1014     float clippedStdev = 0.0;   // self-explanatory
    1015     float oldStanMean = 0.0;    // Temporary variable
    1016     float oldStanStdev = 0.0;   // Temporary variable
     1013    psF32 clippedMean = 0.0;    // self-explanatory
     1014    psF32 clippedStdev = 0.0;   // self-explanatory
     1015    psF32 oldStanMean = 0.0;    // Temporary variable
     1016    psF32 oldStanStdev = 0.0;   // Temporary variable
    10171017    psVector* tmpMask = NULL;   // Temporary vector
    10181018
     
    12401240XXX: Create a 2nd-order polynomial version and solve for X analytically.
    12411241 *****************************************************************************/
    1242 float p_ps1DPolyMedian(psPolynomial1D* myPoly,
    1243                        float rangeLow,
    1244                        float rangeHigh,
    1245                        float getThisValue)
     1242psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
     1243                       psF32 rangeLow,
     1244                       psF32 rangeHigh,
     1245                       psF32 getThisValue)
    12461246{
    12471247    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    12501250    // falls within the range of y-values of the polynomial "myPoly" in the
    12511251    // specified x-range (rangeLow:rangeHigh).
    1252     float fLo = psPolynomial1DEval(
     1252    psF32 fLo = psPolynomial1DEval(
    12531253                    myPoly,
    12541254                    rangeLow
    12551255                );
    1256     float fHi = psPolynomial1DEval(
     1256    psF32 fHi = psPolynomial1DEval(
    12571257                    myPoly,
    12581258                    rangeHigh
     
    12661266
    12671267    psS32 numIterations = 0;
    1268     float midpoint = 0.0;
    1269     float oldMidpoint = 1.0;
    1270     float f = 0.0;
     1268    psF32 midpoint = 0.0;
     1269    psF32 oldMidpoint = 1.0;
     1270    psF32 f = 0.0;
    12711271
    12721272    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
     
    13091309tests to ensure that binNum is within acceptable ranges for both vectors.
    13101310*****************************************************************************/
    1311 float fitQuadraticSearchForYThenReturnX(psVector *xVec,
     1311psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
    13121312                                        psVector *yVec,
    13131313                                        psS32 binNum,
    1314                                         float yVal)
     1314                                        psF32 yVal)
    13151315{
    13161316    PS_VECTOR_CHECK_NULL(xVec, NAN);
     
    13311331    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
    13321332
    1333     float tmpFloat;
     1333    psF32 tmpFloat;
    13341334
    13351335    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
    13361336        // The general case.  We have all three points.
    1337         x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
    1338         x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
    1339         x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
     1337        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
     1338        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
     1339        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
    13401340        y->data.F64[0] = yVec->data.F32[binNum - 1];
    13411341        y->data.F64[1] = yVec->data.F32[binNum];
     
    13461346            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    13471347                    PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
    1348                     (double)yVal,y->data.F64[2],y->data.F64[0]);
     1348                    (psF64)yVal,y->data.F64[2],y->data.F64[0]);
    13491349        }
    13501350        yErr->data.F64[0] = 1.0;
     
    14251425XXX: Check for errors in psLib routines that we call.
    14261426*****************************************************************************/
    1427 int p_psVectorRobustStats(const psVector* restrict myVector,
    1428                           const psVector* restrict errors,
    1429                           const psVector* restrict maskVector,
    1430                           psU32 maskVal,
    1431                           psStats* stats)
     1427psS32 p_psVectorRobustStats(const psVector* restrict myVector,
     1428                            const psVector* restrict errors,
     1429                            const psVector* restrict maskVector,
     1430                            psU32 maskVal,
     1431                            psStats* stats)
    14321432{
    14331433    psHistogram* robustHistogram = NULL;
    14341434    psVector* robustHistogramVector = NULL;
    1435     float binSize = 0.0;        // Size of the histogram bins
     1435    psF32 binSize = 0.0;        // Size of the histogram bins
    14361436    psS32 LQBinNum = -1;          // Bin num for lower quartile
    14371437    psS32 UQBinNum = -1;          // Bin num for upper quartile
     
    14391439    psS32 i = 0;                  // Loop index variable
    14401440    psS32 modeBinNum = 0;
    1441     float modeBinCount = 0.0;
    1442     float dL = 0.0;
     1441    psF32 modeBinCount = 0.0;
     1442    psF32 dL = 0.0;
    14431443    psS32 numBins = 0;
    1444     float myMean = 0.0;
    1445     float myStdev = 0.0;
    1446     float countFloat = 0.0;
    1447     float diff = 0.0;
    1448     float sumSquares = 0.0;
    1449     float sumDiffs = 0.0;
     1444    psF32 myMean = 0.0;
     1445    psF32 myStdev = 0.0;
     1446    psF32 countFloat = 0.0;
     1447    psF32 diff = 0.0;
     1448    psF32 sumSquares = 0.0;
     1449    psF32 sumDiffs = 0.0;
    14501450    psVector* cumulativeRobustSums = NULL;
    1451     float sumRobust = 0.0;
    1452     float sumN50 = 0.0;
    1453     float sumNfit = 0.0;
     1451    psF32 sumRobust = 0.0;
     1452    psF32 sumN50 = 0.0;
     1453    psF32 sumNfit = 0.0;
    14541454    psScalar tmpScalar;
    14551455    tmpScalar.type.type = PS_TYPE_F32;
     
    14601460    // that by 10.0;
    14611461    //XXX: add errors
    1462     int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);
     1462    psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);
    14631463    if (rc != 0) {
    14641464        psError(PS_ERR_UNEXPECTED_NULL,
     
    16141614    for (i=0;i<robustHistogramVector->n;i++) {
    16151615        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1616         ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
     1616        ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    16171617        y->data.F32[i] = robustHistogramVector->data.F32[i];
    16181618    }
     
    17701770    The histogram structure
    17711771 *****************************************************************************/
    1772 psHistogram* psHistogramAlloc(float lower, float upper, psS32 n)
     1772psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n)
    17731773{
    17741774    PS_INT_CHECK_POSITIVE(n, NULL);
     
    17771777    psS32 i = 0;                  // Loop index variable
    17781778    psHistogram* newHist = NULL;        // The new histogram structure
    1779     float binSize = 0.0;        // The histogram bin size
     1779    psF32 binSize = 0.0;        // The histogram bin size
    17801780
    17811781    // Allocate memory for the new histogram structure.  If there are N
     
    17871787
    17881788    // Calculate the bounds for each bin.
    1789     binSize = (upper - lower) / (float)n;
     1789    binSize = (upper - lower) / (psF32)n;
    17901790    // XXX: Is the following necessary? It prevents the max data point
    17911791    // from being in a non-existant bin.
    17921792    binSize += FLT_EPSILON;
    17931793    for (i = 0; i < n + 1; i++) {
    1794         newHist->bounds->data.F32[i] = lower + (binSize * (float)i);
     1794        newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i);
    17951795    }
    17961796
     
    18571857}
    18581858
     1859/*****************************************************************************
     1860UpdateHistogramBins(binNum, out, data, error): This routine is to be used when
     1861updating the histogram in the presence of errors in the input data.  We treat
     1862the data point as a boxcar PDF and update a range of points surrounding the
     1863histogram bin which contains the point.  The width of that boxcar is defined
     1864as 2.35 * error.  Inputs:
     1865    binNum: the bin number of the data point in the histogram
     1866    out: the histogram structure
     1867    data: the data point value
     1868    error: the error in that data point
     1869 
     1870XXX: Must test this.
     1871 *****************************************************************************/
    18591872psS32 UpdateHistogramBins(psS32 binNum,
    18601873                          psHistogram* out,
     
    18631876{
    18641877    PS_PTR_CHECK_NULL(out, -1);
    1865     PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2);
    1866     /*
    1867         psF32 width = 2.35 * error;
    1868         psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum]
    1869         psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0;
    1870      
    1871         if (width <= centerBinWidth) {
    1872             out->nums->data.F32[binNum]+= 1.0;
    1873         } else {
    1874             out->nums->data.F32[binNum]+= centerBinWidth / width;
    1875             // XXX: walk to the left, adding fractional values.
    1876             // XXX: walk to the right, adding fractional values.
    1877      
    1878      
    1879         }
    1880     */
     1878    PS_PTR_CHECK_NULL(out->bounds, -1);
     1879    PS_PTR_CHECK_NULL(out->nums, -1);
     1880    PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2);
     1881    PS_FLOAT_COMPARE(0.0, error, -3);
     1882    PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);
     1883
     1884    psF32 boxcarWidth = 2.35 * error;
     1885    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
     1886                          out->bounds->data.F32[binNum+1]) / 2.0;
     1887    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0);
     1888    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0);
     1889    psS32 bin;
     1890    psS32 boxcarLeftBinNum;
     1891    psS32 boxcarRightBinNum;
     1892
     1893    // Determine the left endpoint of the boxcar for the PDF.
     1894    for (bin=binNum ; bin >= 0 ; bin--) {
     1895        if (out->nums->data.F32[bin] <= boxcarLeft) {
     1896            boxcarLeftBinNum = bin;
     1897            break;
     1898        }
     1899    }
     1900
     1901    // Determine the right endpoint of the boxcar for the PDF.
     1902    for (bin=binNum ; bin < out->nums->n ; bin++) {
     1903        if (out->nums->data.F32[bin] >= boxcarRight) {
     1904            boxcarRightBinNum = bin;
     1905            break;
     1906        }
     1907    }
     1908
     1909    //
     1910    // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
     1911    // bin and return.
     1912    //
     1913    if (boxcarLeftBinNum == boxcarRightBinNum) {
     1914        out->nums->data.F32[binNum]+= 1.0;
     1915        return(0);
     1916    }
     1917
     1918    //
     1919    // If we get here, multiple bins must be updated.  We handle the left
     1920    // endpoint, and right endpoint differently.
     1921    //
     1922    out->nums->data.F32[boxcarLeftBinNum]+=
     1923        (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;
     1924
     1925    //
     1926    // Loop through the center bins, if any.
     1927    //
     1928    for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) {
     1929        out->nums->data.F32[bin]+=
     1930            (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
     1931    }
     1932
     1933    //
     1934    // Handle the right endpoint differently.
     1935    //
     1936    out->nums->data.F32[boxcarRightBinNum]+=
     1937        (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;
     1938
     1939    //
     1940    // Return 0 on success.
     1941    //
    18811942    return(0);
    18821943}
     
    19061967{
    19071968    PS_PTR_CHECK_NULL(out, NULL);
     1969    PS_VECTOR_CHECK_NULL(out->bounds, NULL);
    19081970    PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL);
     1971    PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL);
     1972    PS_VECTOR_CHECK_NULL(out->nums, NULL);
    19091973    PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL);
    19101974    PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL);
     
    19201984
    19211985    psS32 i = 0;                  // Loop index variable
    1922     float binSize = 0.0;          // Histogram bin size
     1986    psF32 binSize = 0.0;          // Histogram bin size
    19231987    psS32 binNum = 0;             // A temporary bin number
    19241988    psS32 numBins = 0;            // The total number of bins
    1925     psS32 tmp = 0;
    19261989    psScalar tmpScalar;
    19271990    tmpScalar.type.type = PS_TYPE_F32;
    1928     psVector* inF32;
    1929     psS32 mustFreeTmp = 1;
     1991    psVector* inF32 = NULL;
     1992    psVector* errorsF32 = NULL;
     1993    psS32 mustFreeVectorIn = 1;
     1994    psS32 mustFreeVectorErrors = 1;
     1995
     1996    // Convert input and errors vectors to F32 if necessary.
    19301997    inF32 = p_psConvertToF32((psVector *) in);
    19311998    if (inF32 == NULL) {
    19321999        inF32 = (psVector *) in;
    1933         mustFreeTmp = 0;
     2000        mustFreeVectorIn = 0;
     2001    }
     2002    errorsF32 = p_psConvertToF32((psVector *) errors);
     2003    if (errorsF32 == NULL) {
     2004        errorsF32 = (psVector *) errors;
     2005        mustFreeVectorErrors = 0;
    19342006    }
    19352007
     
    19502022                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    19512023                    binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
    1952                     if (errors != NULL) {
     2024                    if (errorsF32 != NULL) {
    19532025                        // XXX: Check return codes.
    19542026                        UpdateHistogramBins(binNum, out,
    19552027                                            inF32->data.F32[i],
    1956                                             errors->data.F32[i]);
     2028                                            errorsF32->data.F32[i]);
    19572029                    } else {
    19582030                        // XXX: This if-statement really shouldn't be necessary.
     
    19692041                    // correct bin number requires a bit more work.
    19702042                    tmpScalar.data.F32 = inF32->data.F32[i];
    1971                     tmp = p_psVectorBinDisect(out->bounds, &tmpScalar);
    1972                     if (tmp < 0) {
     2043                    binNum = p_psVectorBinDisect(out->bounds, &tmpScalar);
     2044                    if (binNum < 0) {
    19732045                        psLogMsg(__func__, PS_LOG_WARN,
    19742046                                 "WARNING: psVectorHistogram(): element outside histogram bounds.\n");
    19752047                    } else {
    1976                         if (errors != NULL) {
     2048                        if (errorsF32 != NULL) {
    19772049                            // XXX: Check return codes.
    1978                             UpdateHistogramBins(tmp, out,
     2050                            UpdateHistogramBins(binNum, out,
    19792051                                                inF32->data.F32[i],
    19802052                                                errors->data.F32[i]);
    19812053                        } else {
    1982                             (out->nums->data.F32[tmp])+= 1.0;
     2054                            (out->nums->data.F32[binNum])+= 1.0;
    19832055                        }
    19842056                    }
     
    19882060    }
    19892061
    1990     if (mustFreeTmp == 1) {
     2062    if (mustFreeVectorIn == 1) {
    19912063        psFree(inF32);
     2064    }
     2065    if (mustFreeVectorErrors == 1) {
     2066        psFree(errorsF32);
    19922067    }
    19932068    return (out);
     
    20152090        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20162091        for (i = 0; i < in->n; i++) {
    2017             tmp->data.F32[i] = (float)in->data.S8[i];
     2092            tmp->data.F32[i] = (psF32)in->data.S8[i];
     2093        }
     2094    } else if (in->type.type == PS_TYPE_S16) {
     2095        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2096        for (i = 0; i < in->n; i++) {
     2097            tmp->data.F32[i] = (psF32) in->data.S16[i];
     2098        }
     2099    } else if (in->type.type == PS_TYPE_S32) {
     2100        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2101        for (i = 0; i < in->n; i++) {
     2102            tmp->data.F32[i] = (psF32)in->data.S32[i];
     2103        }
     2104    } else if (in->type.type == PS_TYPE_S64) {
     2105        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2106        for (i = 0; i < in->n; i++) {
     2107            tmp->data.F32[i] = (psF32)in->data.S64[i];
     2108        }
     2109    } else if (in->type.type == PS_TYPE_U8) {
     2110        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2111        for (i = 0; i < in->n; i++) {
     2112            tmp->data.F32[i] = (psF32)in->data.U8[i];
    20182113        }
    20192114    } else if (in->type.type == PS_TYPE_U16) {
    20202115        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20212116        for (i = 0; i < in->n; i++) {
    2022             tmp->data.F32[i] = (float)in->data.U16[i];
    2023         }
    2024     } else if (in->type.type == PS_TYPE_U8) {
     2117            tmp->data.F32[i] = (psF32)in->data.U16[i];
     2118        }
     2119    } else if (in->type.type == PS_TYPE_U32) {
    20252120        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20262121        for (i = 0; i < in->n; i++) {
    2027             tmp->data.F32[i] = (float)in->data.U8[i];
     2122            tmp->data.F32[i] = (psF32)in->data.U32[i];
     2123        }
     2124    } else if (in->type.type == PS_TYPE_U64) {
     2125        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2126        for (i = 0; i < in->n; i++) {
     2127            tmp->data.F32[i] = (psF32)in->data.U64[i];
    20282128        }
    20292129    } else if (in->type.type == PS_TYPE_F64) {
    20302130        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20312131        for (i = 0; i < in->n; i++) {
    2032             tmp->data.F32[i] = (float)in->data.F64[i];
     2132            tmp->data.F32[i] = (psF32)in->data.F64[i];
    20332133        }
    20342134    } else if (in->type.type == PS_TYPE_F32) {
    20352135        // do nothing
    20362136    } else {
    2037         char* strType;
    2038         PS_TYPE_NAME(strType,in->type.type);
     2137        psS8* strType;
     2138        PS_TYPE_NAME(strType, in->type.type);
    20392139        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    20402140                PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED,
     
    20742174    }
    20752175
    2076     psVector* inF32;
    2077     psVector* errorsF32;
     2176    psVector* inF32 = NULL;
     2177    psVector* errorsF32 = NULL;
    20782178    psS32 mustFreeVectorIn = 1;
    20792179    psS32 mustFreeVectorErrors = 1;
     
    21782278    return (stats);
    21792279}
     2280
  • trunk/psLib/src/dataManip/psStats.h

    r2778 r2788  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-12-21 20:42:07 $
     12 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-12-22 05:09:32 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5656typedef struct
    5757{
    58     double sampleMean;          ///< formal mean of sample
    59     double sampleMedian;        ///< formal median of sample
    60     double sampleStdev;         ///< standard deviation of sample
    61     double sampleUQ;            ///< upper quartile of sample
    62     double sampleLQ;            ///< lower quartile of sample
    63     double robustMean;          ///< robust mean of array
    64     double robustMedian;        ///< robust median of array
    65     double robustMode;          ///< Robust mode of array
    66     double robustStdev;         ///< robust standard deviation of array
    67     double robustUQ;            ///< robust upper quartile
    68     double robustLQ;            ///< robust lower quartile
     58    psF64 sampleMean;          ///< formal mean of sample
     59    psF64 sampleMedian;        ///< formal median of sample
     60    psF64 sampleStdev;         ///< standard deviation of sample
     61    psF64 sampleUQ;            ///< upper quartile of sample
     62    psF64 sampleLQ;            ///< lower quartile of sample
     63    psF64 robustMean;          ///< robust mean of array
     64    psF64 robustMedian;        ///< robust median of array
     65    psF64 robustMode;          ///< Robust mode of array
     66    psF64 robustStdev;         ///< robust standard deviation of array
     67    psF64 robustUQ;            ///< robust upper quartile
     68    psF64 robustLQ;            ///< robust lower quartile
    6969    psS32 robustN50;              ///<
    7070    psS32 robustNfit;             ///<
    71     double clippedMean;         ///< Nsigma clipped mean
    72     double clippedStdev;        ///< standard deviation after clipping
     71    psF64 clippedMean;         ///< Nsigma clipped mean
     72    psF64 clippedStdev;        ///< standard deviation after clipping
    7373    psS32 clippedNvalues;         ///< ???
    74     double clipSigma;           ///< Nsigma used for clipping; user input
     74    psF64 clipSigma;           ///< Nsigma used for clipping; user input
    7575    psS32 clipIter;               ///< Number of clipping iterations; user input
    76     double min;                 ///< minimum data value in array
    77     double max;                 ///< maximum data value in array
    78     double binsize;             ///<
     76    psF64 min;                 ///< minimum data value in array
     77    psF64 max;                 ///< maximum data value in array
     78    psF64 binsize;             ///<
    7979    psStatsOptions options;     ///< bitmask of calculated values
    8080}
     
    132132 */
    133133psHistogram* psHistogramAlloc(
    134     float lower,                       ///< Lower limit for the bins
    135     float upper,                       ///< Upper limit for the bins
     134    psF32 lower,                       ///< Lower limit for the bins
     135    psF32 upper,                       ///< Upper limit for the bins
    136136    psS32 n                              ///< Number of bins
    137137);
     
    172172    ///< the statistic struct to operate on
    173173
    174     double *value
     174    psF64 *value
    175175    ///< if return is true, this is set to the specified statistic value by stats->options
    176176);
  • trunk/psLib/src/math/psConstants.h

    r2778 r2788  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-12-21 20:42:07 $
     8 *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-12-22 05:09:32 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     12 *
     13 *  XXX: Add parenthesis around all arguments so that these macros can be
     14 *       called with complex expressions.
     15 *
     16 *
     17 *
    1218 */
    1319
  • trunk/psLib/src/math/psMinimize.c

    r2741 r2788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-17 00:18:31 $
     11 *  @version $Revision: 1.96 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-22 05:09:32 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 *
    1616 *  XXX: must follow coding name standards on local functions.
    17  *
    18  *  XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via
    19  *       LMM) is not addressed here.  We are waiting for subsequent SDRs
    20  *       which will redefine the LMM functions.
    2117 *
    2218 */
     
    6460XXX: Use a static vector.
    6561 *****************************************************************************/
    66 void psBuildSums1D(double x,
     62void psBuildSums1D(psF64 x,
    6763                   psS32 polyOrder,
    6864                   psVector* sums)
    6965{
    7066    psS32 i = 0;
    71     double xSum = 0.0;
     67    psF64 xSum = 0.0;
    7268
    7369    if (sums == NULL) {
     
    10096XXX: do an F64 version?
    10197 *****************************************************************************/
    102 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
     98psF32 *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
    10399                             const psVector* restrict y)        ///< Coordinates
    104100{
     
    108104    psS32 i;
    109105    psS32 k;
    110     float sig;
    111     float p;
     106    psF32 sig;
     107    psF32 p;
    112108    psS32 n = y->n;
    113     float *u = (float *) psAlloc(n * sizeof(float));
    114     float *derivs2 = (float *) psAlloc(n * sizeof(float));
    115     float *X = (float *) & (x->data.F32[0]);
    116     float *Y = (float *) & (y->data.F32[0]);
    117     float qn;
     109    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
     110    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
     111    psF32 *X = (psF32 *) & (x->data.F32[0]);
     112    psF32 *Y = (psF32 *) & (y->data.F32[0]);
     113    psF32 qn;
    118114
    119115    // XXX: The second derivatives at the endpoints, undefined in the SDR,
     
    168164 *****************************************************************************/
    169165/*
    170 float p_psNRSpline1DEval(psSpline1D *spline,
     166psF32 p_psNRSpline1DEval(psSpline1D *spline,
    171167                         const psVector* restrict x,
    172168                         const psVector* restrict y,
    173                          float X)
     169                         psF32 X)
    174170{
    175171    PS_PTR_CHECK_NULL(spline, NAN);
     
    185181    psS32 klo;
    186182    psS32 khi;
    187     float H;
    188     float A;
    189     float B;
    190     float C;
    191     float D;
    192     float Y;
     183    psF32 H;
     184    psF32 A;
     185    psF32 B;
     186    psF32 C;
     187    psF32 D;
     188    psF32 Y;
    193189 
    194190    n = spline->n;
     
    258254            "---- psVectorFitSpline1D() begin ----\n");
    259255    psS32 numSplines = (y->n)-1;
    260     float tmp;
    261     float H;
     256    psF32 tmp;
     257    psF32 H;
    262258    psS32 i;
    263     float slope;
     259    psF32 slope;
    264260    psVector *x32 = NULL;
    265261    psVector *y32 = NULL;
     
    421417    psTrace(".psLib.dataManip.psMinimize", 4,
    422418            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
    423     float x;
    424     int i;
    425     float mean = params->data.F32[0];
    426     float stdev = params->data.F32[1];
     419    psF32 x;
     420    psS32 i;
     421    psF32 mean = params->data.F32[0];
     422    psF32 stdev = params->data.F32[1];
    427423    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    428424
     
    443439    for (i=0;i<coords->n;i++) {
    444440        x = ((psVector *) (coords->data[i]))->data.F32[0];
    445         float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
     441        psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    446442        deriv->data.F32[i][0] = tmp / (stdev * stdev);
    447443        tmp = (x - mean) * (x - mean) *
     
    491487    PS_PTR_CHECK_NULL(params, NULL);
    492488
    493     double normalization = params->data.F32[0];
    494     double x0 = params->data.F32[1];
    495     double y0 = params->data.F32[2];
    496     double sigmaX = params->data.F32[3];
    497     double sigmaY = params->data.F32[4];
    498     double theta = params->data.F32[5];
     489    psF64 normalization = params->data.F32[0];
     490    psF64 x0 = params->data.F32[1];
     491    psF64 y0 = params->data.F32[2];
     492    psF64 sigmaX = params->data.F32[3];
     493    psF64 sigmaY = params->data.F32[4];
     494    psF64 theta = params->data.F32[5];
    499495    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    500496
     
    502498        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
    503499    } else {
    504         // XXX: Check size of derivative
     500        PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL);
    505501    }
    506502
     
    508504            "---- psMinimizeLMChi2Gauss2D() begin ----\n");
    509505
    510     for (int i=0;i<coords->n;i++) {
    511         double x = ((psVector *) coords->data[i])->data.F32[0];
    512         double y = ((psVector *) coords->data[i])->data.F32[0];
    513 
    514         double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
    515         double v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
    516 
    517         double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
    518                                              v*v/(2.0 * sigmaY * sigmaY)))/
    519                       (2.0 * M_PI * sigmaX * sigmaY);
     506    for (psS32 i=0;i<coords->n;i++) {
     507        psF64 x = ((psVector *) coords->data[i])->data.F32[0];
     508        psF64 y = ((psVector *) coords->data[i])->data.F32[0];
     509
     510        psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta);
     511        psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta);
     512
     513        psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +
     514                                            v*v/(2.0 * sigmaY * sigmaY)))/
     515                     (2.0 * M_PI * sigmaX * sigmaY);
    520516        out->data.F32[i] = flux;
    521517
     
    584580        PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL);
    585581    }
    586 
    587     // XXX: Generate code that modifies covar matrix if not-NULL.
    588582    if (covar != NULL) {
    589583        PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL);
     
    610604    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    611605    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    612 
    613     //    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
    614     //    for (i=0;i<numData;i++) {
    615     //        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
    616     //    }
    617606    psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32);
    618 
    619     psVector *currValueVec;
    620     psVector *newValueVec;
    621 
    622     float currChi2 = 0.0;
    623     float newChi2 = 0.0;
    624     float lamda = 0.00005;
     607    psVector *currValueVec = NULL;
     608    psVector *newValueVec = NULL;
     609    psF32 currChi2 = 0.0;
     610    psF32 newChi2 = 0.0;
     611    psF32 lamda = 0.00005;
    625612    lamda = 0.05;
    626613
     
    711698                if (j == k) {
    712699                    A->data.F64[j][k] =
    713                         (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
     700                        (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]);
    714701                } else {
    715                     A->data.F64[j][k] = (double) alpha->data.F32[j][k];
     702                    A->data.F64[j][k] = (psF64) alpha->data.F32[j][k];
    716703                }
    717704            }
     
    741728            } else {
    742729                newParams->data.F32[i] = params->data.F32[i] -
    743                                          (float) paramDeltasF64->data.F64[i];
     730                                         (psF32) paramDeltasF64->data.F64[i];
    744731            }
    745732        }
     
    777764            // We already masked params.
    778765            for (i=0;i<numParams;i++) {
    779                 params->data.F32[i] = (float) newParams->data.F32[i];
     766                params->data.F32[i] = (psF32) newParams->data.F32[i];
    780767            }
    781768            lamda*= 0.1;
     
    826813    psS32 k;
    827814    psS32 n = x->n;
    828     double fac;
    829     double sum;
     815    psF64 fac;
     816    psF64 sum;
    830817    PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64);
    831818    psScalar *fScalar;
     
    836823    // variable declarations.  I retain them here to maintain coherence
    837824    // with the NR code.
    838     double min = -1.0;
    839     double max = 1.0;
    840     double bma = 0.5 * (max-min);  // 1
    841     double bpa = 0.5 * (max+min);  // 0
     825    psF64 min = -1.0;
     826    psF64 max = 1.0;
     827    psF64 bma = 0.5 * (max-min);  // 1
     828    psF64 bpa = 0.5 * (max+min);  // 0
    842829
    843830    // In this loop, we first calculate the values of X for which the
     
    851838    for (psS32 i=0;i<n;i++) {
    852839        // NR 5.8.4
    853         double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
    854         double X = (Y + bma + bpa) - 1.0;
     840        psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n));
     841        psF64 X = (Y + bma + bpa) - 1.0;
    855842        tmpScalar.data.F64 = X;
    856843
     
    873860    // coefficients of the Chebyshev polynomial: NR 5.8.7.
    874861
    875     fac = 2.0/((float) n);
     862    fac = 2.0/((psF32) n);
    876863    // XXX: is this loop bound correct?
    877864    for (j=0;j<myPoly->n;j++) {
     
    879866        for (k=0;k<n;k++) {
    880867            sum+= f->data.F64[k] *
    881                   cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
     868                  cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n));
    882869        }
    883870
     
    10751062 *****************************************************************************/
    10761063psMinimization *psMinimizationAlloc(psS32 maxIter,
    1077                                     float tol)
     1064                                    psF32 tol)
    10781065{
    10791066    PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL);
     
    10921079// LINE to it.  We assume BASEMASK is non-null.
    10931080#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
    1094 for (int i=0;i<BASE->n;i++) { \
     1081for (psS32 i=0;i<BASE->n;i++) { \
    10951082    if (BASEMASK->data.U8[i] == 0) { \
    10961083        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
     
    11021089#define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \
    11031090BOOL_VAR = true; \
    1104 for (int i=0;i<IN->n;i++) { \
     1091for (psS32 i=0;i<IN->n;i++) { \
    11051092    if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \
    11061093        BOOL_VAR = false; \
     
    11111098#define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \
    11121099BOOL_VAR = true; \
    1113 for (int i=0;i<IN->n;i++) { \
     1100for (psS32 i=0;i<IN->n;i++) { \
    11141101    if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \
    11151102        BOOL_VAR = false; \
     
    11521139                               psMinimizePowellFunc func)
    11531140{
    1154     float a = 0.0;
    1155     float b = 0.0;
    1156     float c = 0.0;
    1157     float fa = 0.0;
    1158     float fb = 0.0;
    1159     float fc = 0.0;
     1141    psF32 a = 0.0;
     1142    psF32 b = 0.0;
     1143    psF32 c = 0.0;
     1144    psF32 fa = 0.0;
     1145    psF32 fb = 0.0;
     1146    psF32 fc = 0.0;
    11601147    psS32 iter = 100;
    1161     float aDir = 0.0;
    1162     float cDir = 0.0;
    1163     float new_aDir = 0.0;
    1164     float new_cDir = 0.0;
     1148    psF32 aDir = 0.0;
     1149    psF32 cDir = 0.0;
     1150    psF32 new_aDir = 0.0;
     1151    psF32 new_cDir = 0.0;
    11651152    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1166     float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
     1153    psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
    11671154    psVector *tmp = NULL;
    11681155    psBool boolLineIsNull = true;
     
    13011288                                psMinimizePowellFunc func)
    13021289{
    1303     float a = 0.0;
    1304     float b = 0.0;
    1305     float c = 0.0;
    1306     float fa = 0.0;
    1307     float fb = 0.0;
    1308     float fc = 0.0;
     1290    psF32 a = 0.0;
     1291    psF32 b = 0.0;
     1292    psF32 c = 0.0;
     1293    psF32 fa = 0.0;
     1294    psF32 fb = 0.0;
     1295    psF32 fc = 0.0;
    13091296    psS32 iter = 0;
    13101297    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
    13111298    psBool boolLineIsNull = true;
    1312     float prevMin = 0.0;
    1313     int countMin = 0;
     1299    psF32 prevMin = 0.0;
     1300    psS32 countMin = 0;
    13141301
    13151302    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
     
    14031390 *****************************************************************************/
    14041391#define PS_LINEMIN_MAX_ITERATIONS 30
    1405 float p_psLineMin(psMinimization *min,
     1392psF32 p_psLineMin(psMinimization *min,
    14061393                  psVector *params,
    14071394                  psVector *line,
     
    14231410    PS_PTR_CHECK_NULL(func, NAN);
    14241411    psVector *bracket;
    1425     float a = 0.0;
    1426     float b = 0.0;
    1427     float c = 0.0;
    1428     float n = 0.0;
    1429     float fa = 0.0;
    1430     float fb = 0.0;
    1431     float fc = 0.0;
    1432     float fn = 0.0;
    1433     float mul = 0.0;
     1412    psF32 a = 0.0;
     1413    psF32 b = 0.0;
     1414    psF32 c = 0.0;
     1415    psF32 n = 0.0;
     1416    psF32 fa = 0.0;
     1417    psF32 fb = 0.0;
     1418    psF32 fc = 0.0;
     1419    psF32 fn = 0.0;
     1420    psF32 mul = 0.0;
    14341421    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
    14351422    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
     
    15441531This routine must minimize a possibly multi-dimensional function.  The
    15451532function to be minimized "func" is:
    1546     float func(psVector *params, psArray *coords)
     1533    psF32 func(psVector *params, psArray *coords)
    15471534The "params" are the parameters of the function which are varied.  The data
    15481535points at which the function is varied are in the argument "coords" which is
     
    15781565    psVector *myParamMask = NULL;
    15791566    psMinimization dummyMin;
    1580     float mul = 0.0;
    1581     float baseFuncVal = 0.0;
    1582     float currFuncVal = 0.0;
     1567    psF32 mul = 0.0;
     1568    psF32 baseFuncVal = 0.0;
     1569    psF32 currFuncVal = 0.0;
    15831570    psS32 biggestIter = 0;
    1584     float biggestDiff = 0.0;
    1585     int iterationNumber = 0;
     1571    psF32 biggestDiff = 0.0;
     1572    psS32 iterationNumber = 0;
    15861573
    15871574    psTrace(".psLib.dataManip.psMinimizePowell", 4,
     
    17111698            }
    17121699        }
    1713         float fqp = func(pQP, coords);
    1714         float term1 = (baseFuncVal - currFuncVal) - biggestDiff;
     1700        psF32 fqp = func(pQP, coords);
     1701        psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff;
    17151702        term1*= term1;
    17161703        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
    1717         float term2 = baseFuncVal - fqp;
     1704        psF32 term2 = baseFuncVal - fqp;
    17181705        term2*= term2 * biggestDiff;
    17191706        if (term1 < term2) {
     
    17601747    PS_PTR_CHECK_NULL(params, NULL);
    17611748
    1762     float x;
    1763     int i;
    1764     float mean = params->data.F32[0];
    1765     float stdev = params->data.F32[1];
     1749    psF32 x;
     1750    psS32 i;
     1751    psF32 mean = params->data.F32[0];
     1752    psF32 stdev = params->data.F32[1];
    17661753    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
    17671754
     
    17831770XXX: This is F32 only
    17841771 *****************************************************************************/
    1785 float myPowellChi2Func(const psVector *params,
     1772psF32 myPowellChi2Func(const psVector *params,
    17861773                       const psArray *coords)
    17871774{
     
    17941781    PS_PTR_CHECK_NULL(coords, NAN);
    17951782
    1796     float chi2 = 0.0;
    1797     float d;
     1783    psF32 chi2 = 0.0;
     1784    psF32 d;
    17981785    psS32 i;
    17991786    psVector *tmp;
  • trunk/psLib/src/math/psMinimize.h

    r2250 r2788  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-11-01 23:57:08 $
     10 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-12-22 05:09:32 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4040{
    4141    psS32 maxIter;                       ///< Convergence limit
    42     float tol;                         ///< Error Tolerance
    43     float value;                       ///< Value of function at minimum
     42    psF32 tol;                         ///< Error Tolerance
     43    psF32 value;                       ///< Value of function at minimum
    4444    psS32 iter;                          ///< Number of iterations required
    45     float lastDelta;                   ///< The last difference for the fit
     45    psF32 lastDelta;                   ///< The last difference for the fit
    4646}
    4747psMinimization;
    4848
    4949psMinimization *psMinimizationAlloc(psS32 maxIter,
    50                                     float tol);
     50                                    psF32 tol);
    5151
    5252/** Derive a polynomial fit.
     
    9191
    9292typedef
    93 float (*psMinimizePowellFunc)(const psVector *params,
     93psF32 (*psMinimizePowellFunc)(const psVector *params,
    9494                              const psArray *coords);
    9595
  • trunk/psLib/src/math/psPolynomial.c

    r2778 r2788  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-12-21 20:42:07 $
     9 *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-12-22 05:09:32 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5555static void dPolynomial4DFree(psDPolynomial4D* myPoly);
    5656static void spline1DFree(psSpline1D *tmpSpline);
    57 static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x);
     57static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
    5858static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
    5959
     
    277277 
    278278 *****************************************************************************/
    279 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     279static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    280280{
    281281    psS32 loop_x = 0;
    282     float polySum = 0.0;
    283     float xSum = 1.0;
     282    psF32 polySum = 0.0;
     283    psF32 xSum = 1.0;
    284284
    285285    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
     
    307307// XXX: How does the mask vector effect Crenshaw's formula?
    308308// XXX: We assume that x is scaled between -1.0 and 1.0;
    309 static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     309static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    310310{
    311311    //    PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0);
     
    314314    psS32 n;
    315315    psS32 i;
    316     float tmp;
     316    psF32 tmp;
    317317
    318318    n = myPoly->n;
     
    337337    psS32 n;
    338338    psS32 i;
    339     float tmp;
     339    psF32 tmp;
    340340    psPolynomial1D **chebPolys = NULL;
    341341
     
    355355}
    356356
    357 static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     357static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    358358{
    359359    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    361361    psS32 loop_x = 0;
    362362    psS32 loop_y = 0;
    363     float polySum = 0.0;
    364     float xSum = 1.0;
    365     float ySum = 1.0;
     363    psF32 polySum = 0.0;
     364    psF32 xSum = 1.0;
     365    psF32 ySum = 1.0;
    366366
    367367    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    379379}
    380380
    381 static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     381static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    382382{
    383383    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    386386    psS32 loop_y = 0;
    387387    psS32 i = 0;
    388     float polySum = 0.0;
     388    psF32 polySum = 0.0;
    389389    psPolynomial1D* *chebPolys = NULL;
    390390    psS32 maxChebyPoly = 0;
     
    414414}
    415415
    416 static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     416static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    417417{
    418418    psS32 loop_x = 0;
    419419    psS32 loop_y = 0;
    420420    psS32 loop_z = 0;
    421     float polySum = 0.0;
    422     float xSum = 1.0;
    423     float ySum = 1.0;
    424     float zSum = 1.0;
     421    psF32 polySum = 0.0;
     422    psF32 xSum = 1.0;
     423    psF32 ySum = 1.0;
     424    psF32 zSum = 1.0;
    425425
    426426    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    442442}
    443443
    444 static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     444static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    445445{
    446446    psS32 loop_x = 0;
     
    448448    psS32 loop_z = 0;
    449449    psS32 i = 0;
    450     float polySum = 0.0;
     450    psF32 polySum = 0.0;
    451451    psPolynomial1D* *chebPolys = NULL;
    452452    psS32 maxChebyPoly = 0;
     
    483483}
    484484
    485 static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     485static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    486486{
    487487    psS32 loop_w = 0;
     
    489489    psS32 loop_y = 0;
    490490    psS32 loop_z = 0;
    491     float polySum = 0.0;
    492     float wSum = 1.0;
    493     float xSum = 1.0;
    494     float ySum = 1.0;
    495     float zSum = 1.0;
     491    psF32 polySum = 0.0;
     492    psF32 wSum = 1.0;
     493    psF32 xSum = 1.0;
     494    psF32 ySum = 1.0;
     495    psF32 zSum = 1.0;
    496496
    497497    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    517517}
    518518
    519 static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     519static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    520520{
    521521    psS32 loop_w = 0;
     
    524524    psS32 loop_z = 0;
    525525    psS32 i = 0;
    526     float polySum = 0.0;
     526    psF32 polySum = 0.0;
    527527    psPolynomial1D* *chebPolys = NULL;
    528528    psS32 maxChebyPoly = 0;
     
    568568    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    569569 *****************************************************************************/
    570 static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     570static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    571571{
    572572    psS32 loop_x = 0;
    573     double polySum = 0.0;
    574     double xSum = 1.0;
     573    psF64 polySum = 0.0;
     574    psF64 xSum = 1.0;
    575575
    576576    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     
    586586// XXX: You can do this without having to psAlloc() vector d.
    587587// XXX: How does the mask vector effect Crenshaw's formula?
    588 static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     588static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    589589{
    590590    psVector *d;
    591591    psS32 n;
    592592    psS32 i;
    593     double tmp;
     593    psF64 tmp;
    594594
    595595    n = myPoly->n;
     
    611611}
    612612
    613 static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     613static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    614614{
    615615    psS32 loop_x = 0;
    616616    psS32 loop_y = 0;
    617     double polySum = 0.0;
    618     double xSum = 1.0;
    619     double ySum = 1.0;
     617    psF64 polySum = 0.0;
     618    psF64 xSum = 1.0;
     619    psF64 ySum = 1.0;
    620620
    621621    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    633633}
    634634
    635 static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     635static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    636636{
    637637    psS32 loop_x = 0;
    638638    psS32 loop_y = 0;
    639639    psS32 i = 0;
    640     double polySum = 0.0;
     640    psF64 polySum = 0.0;
    641641    psPolynomial1D* *chebPolys = NULL;
    642642    psS32 maxChebyPoly = 0;
     
    667667}
    668668
    669 static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     669static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    670670{
    671671    psS32 loop_x = 0;
    672672    psS32 loop_y = 0;
    673673    psS32 loop_z = 0;
    674     double polySum = 0.0;
    675     double xSum = 1.0;
    676     double ySum = 1.0;
    677     double zSum = 1.0;
     674    psF64 polySum = 0.0;
     675    psF64 xSum = 1.0;
     676    psF64 ySum = 1.0;
     677    psF64 zSum = 1.0;
    678678
    679679    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    695695}
    696696
    697 static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     697static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    698698{
    699699    psS32 loop_x = 0;
     
    701701    psS32 loop_z = 0;
    702702    psS32 i = 0;
    703     double polySum = 0.0;
     703    psF64 polySum = 0.0;
    704704    psPolynomial1D* *chebPolys = NULL;
    705705    psS32 maxChebyPoly = 0;
     
    736736}
    737737
    738 static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     738static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    739739{
    740740    psS32 loop_w = 0;
     
    742742    psS32 loop_y = 0;
    743743    psS32 loop_z = 0;
    744     double polySum = 0.0;
    745     double wSum = 1.0;
    746     double xSum = 1.0;
    747     double ySum = 1.0;
    748     double zSum = 1.0;
     744    psF64 polySum = 0.0;
     745    psF64 wSum = 1.0;
     746    psF64 xSum = 1.0;
     747    psF64 ySum = 1.0;
     748    psF64 zSum = 1.0;
    749749
    750750    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    770770}
    771771
    772 static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     772static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    773773{
    774774    psS32 loop_w = 0;
     
    777777    psS32 loop_z = 0;
    778778    psS32 i = 0;
    779     double polySum = 0.0;
     779    psF64 polySum = 0.0;
    780780    psPolynomial1D* *chebPolys = NULL;
    781781    psS32 maxChebyPoly = 0;
     
    828828 *****************************************************************************/
    829829#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
    830 static float fullInterpolate1D##TYPE(ps##TYPE *domain, \
     830static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
    831831                                     ps##TYPE *range, \
    832832                                     psS32 n, \
     
    896896LaGrange interpolation.
    897897 *****************************************************************************/
    898 static float interpolate1DF32(float *domain,
    899                               float *range,
     898static psF32 interpolate1DF32(psF32 *domain,
     899                              psF32 *range,
    900900                              psS32 n,
    901901                              psS32 order,
    902                               float x)
     902                              psF32 x)
    903903{
    904904    psS32 binNum;
     
    941941    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    942942 *****************************************************************************/
    943 float psGaussian(float x, float mean, float sigma, psBool normal)
    944 {
    945     float tmp = 1.0;
     943psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal)
     944{
     945    psF32 tmp = 1.0;
    946946
    947947    psTrace(".psLib.dataManip.psFunctions.psGaussian", 4,
     
    967967XXX: There is no way to seed the random generator.
    968968 *****************************************************************************/
    969 psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts)
     969psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
    970970{
    971971    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     
    10071007    newPoly->type = type;
    10081008    newPoly->n = n;
    1009     newPoly->coeff = (float *)psAlloc(n * sizeof(float));
    1010     newPoly->coeffErr = (float *)psAlloc(n * sizeof(float));
    1011     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1009    newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32));
     1010    newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32));
     1011    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    10121012    for (i = 0; i < n; i++) {
    10131013        newPoly->coeff[i] = 0.0;
     
    10361036    newPoly->nY = nY;
    10371037
    1038     newPoly->coeff = (float **)psAlloc(nX * sizeof(float *));
    1039     newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *));
    1040     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1038    newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1039    newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1040    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    10411041    for (x = 0; x < nX; x++) {
    1042         newPoly->coeff[x] = (float *)psAlloc(nY * sizeof(float));
    1043         newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float));
    1044         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1042        newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1043        newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1044        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    10451045    }
    10461046    for (x = 0; x < nX; x++) {
     
    10751075    newPoly->nZ = nZ;
    10761076
    1077     newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **));
    1078     newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **));
    1079     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1077    newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1078    newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1079    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    10801080    for (x = 0; x < nX; x++) {
    1081         newPoly->coeff[x] = (float **)psAlloc(nY * sizeof(float *));
    1082         newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *));
    1083         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1081        newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1082        newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1083        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    10841084        for (y = 0; y < nY; y++) {
    1085             newPoly->coeff[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1086             newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1087             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1085            newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1086            newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1087            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    10881088        }
    10891089    }
     
    11241124    newPoly->nZ = nZ;
    11251125
    1126     newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***));
    1127     newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***));
    1128     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1126    newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1127    newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1128    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    11291129    for (w = 0; w < nW; w++) {
    1130         newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **));
    1131         newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **));
    1132         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1130        newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1131        newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1132        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    11331133        for (x = 0; x < nX; x++) {
    1134             newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1135             newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1136             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1134            newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1135            newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1136            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    11371137            for (y = 0; y < nY; y++) {
    1138                 newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1139                 newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1140                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1138                newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1139                newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1140                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    11411141            }
    11421142        }
     
    11571157}
    11581158
    1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, float x)
     1159psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x)
    11601160{
    11611161    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    11951195}
    11961196
    1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y)
     1197psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y)
    11981198{
    11991199    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    12501250}
    12511251
    1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z)
     1252psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z)
    12531253{
    12541254    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    13161316}
    13171317
    1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z)
     1318psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z)
    13191319{
    13201320    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    14061406    newPoly->type = type;
    14071407    newPoly->n = n;
    1408     newPoly->coeff = (double *)psAlloc(n * sizeof(double));
    1409     newPoly->coeffErr = (double *)psAlloc(n * sizeof(double));
    1410     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1408    newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64));
     1409    newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64));
     1410    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    14111411    for (i = 0; i < n; i++) {
    14121412        newPoly->coeff[i] = 0.0;
     
    14351435    newPoly->nY = nY;
    14361436
    1437     newPoly->coeff = (double **)psAlloc(nX * sizeof(double *));
    1438     newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *));
    1439     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1437    newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1438    newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1439    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    14401440    for (x = 0; x < nX; x++) {
    1441         newPoly->coeff[x] = (double *)psAlloc(nY * sizeof(double));
    1442         newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double));
    1443         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1441        newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1442        newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1443        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    14441444    }
    14451445    for (x = 0; x < nX; x++) {
     
    14741474    newPoly->nZ = nZ;
    14751475
    1476     newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **));
    1477     newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **));
    1478     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1476    newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1477    newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1478    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    14791479    for (x = 0; x < nX; x++) {
    1480         newPoly->coeff[x] = (double **)psAlloc(nY * sizeof(double *));
    1481         newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *));
    1482         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1480        newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1481        newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1482        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    14831483        for (y = 0; y < nY; y++) {
    1484             newPoly->coeff[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1485             newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1486             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1484            newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1485            newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1486            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    14871487        }
    14881488    }
     
    15231523    newPoly->nZ = nZ;
    15241524
    1525     newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***));
    1526     newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***));
    1527     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1525    newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1526    newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1527    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    15281528    for (w = 0; w < nW; w++) {
    1529         newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **));
    1530         newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **));
    1531         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1529        newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1530        newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1531        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    15321532        for (x = 0; x < nX; x++) {
    1533             newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1534             newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1535             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1533            newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1534            newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1535            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    15361536            for (y = 0; y < nY; y++) {
    1537                 newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1538                 newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1539                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1537                newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1538                newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1539                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    15401540            }
    15411541        }
     
    15571557
    15581558
    1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x)
     1559psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x)
    15601560{
    15611561    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    15981598
    15991599
    1600 double psDPolynomial2DEval(const psDPolynomial2D* myPoly,
    1601                            double x,
    1602                            double y)
     1600psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly,
     1601                          psF64 x,
     1602                          psF64 y)
    16031603{
    16041604    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    16551655
    16561656
    1657 double psDPolynomial3DEval(const psDPolynomial3D* myPoly,
    1658                            double x,
    1659                            double y,
    1660                            double z)
     1657psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly,
     1658                          psF64 x,
     1659                          psF64 y,
     1660                          psF64 z)
    16611661{
    16621662    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    17241724}
    17251725
    1726 double psDPolynomial4DEval(const psDPolynomial4D* myPoly,
    1727                            double w,
    1728                            double x,
    1729                            double y,
    1730                            double z)
     1726psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly,
     1727                          psF64 w,
     1728                          psF64 x,
     1729                          psF64 y,
     1730                          psF64 z)
    17311731{
    17321732    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    18101810//    psS32 n;
    18111811//    psPolynomial1D **spline;
    1812 //    float *p_psDeriv2;
    1813 //    float *domains;
     1812//    psF32 *p_psDeriv2;
     1813//    psF32 *domains;
    18141814//} psSpline1D;
    18151815
     
    18221822psSpline1D *psSpline1DAlloc(psS32 numSplines,
    18231823                            psS32 order,
    1824                             float min,
    1825                             float max)
     1824                            psF32 min,
     1825                            psF32 max)
    18261826{
    18271827    PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL);
     
    18311831    psSpline1D *tmp = NULL;
    18321832    psS32 i;
    1833     float tmpDomain;
    1834     float width;
     1833    psF32 tmpDomain;
     1834    psF32 width;
    18351835
    18361836    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
     
    18451845    tmp->p_psDeriv2 = NULL;
    18461846
    1847     tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    1848     width = (max - min) / ((float) numSplines);
     1847    tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32));
     1848    width = (max - min) / ((psF32) numSplines);
    18491849
    18501850    (tmp->domains)[0] = min;
     
    18841884    }
    18851885
    1886     tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float));
     1886    tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));
    18871887
    18881888    for (i=0;i<bounds->n;i++) {
     
    19771977{
    19781978    PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3);
    1979     char* strType;
     1979    psS8* strType;
    19801980
    19811981    switch (x->type.type) {
     
    20742074        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
    20752075
    2076         psScalar *tmpScalar = psScalarAlloc((double)
     2076        psScalar *tmpScalar = psScalarAlloc((psF64)
    20772077                                            interpolate1DF32(domain32->data.F32,
    20782078                                                             range32->data.F32,
    20792079                                                             domain32->n,
    20802080                                                             order,
    2081                                                              (float) x->data.F64), PS_TYPE_F64);
     2081                                                             (psF32) x->data.F64), PS_TYPE_F64);
    20822082        psFree(range32);
    20832083        psFree(domain32);
     
    20892089
    20902090    } else {
    2091         char* strType;
     2091        psS8* strType;
    20922092        PS_TYPE_NAME(strType,x->type.type);
    20932093        psError(PS_ERR_BAD_PARAMETER_TYPE,
     
    21152115     the spline fit functions require F32 and F64.
    21162116 *****************************************************************************/
    2117 float psSpline1DEval(
    2118     float x,
     2117psF32 psSpline1DEval(
     2118    psF32 x,
    21192119    const psSpline1D *spline
    21202120)
     
    21722172        for (i=0;i<x->n;i++) {
    21732173            tmpVector->data.F32[i] = psSpline1DEval(
    2174                                          (float) x->data.F64[i],
     2174                                         (psF32) x->data.F64[i],
    21752175                                         spline
    21762176                                     );
    21772177        }
    21782178    } else {
    2179         char* strType;
     2179        psS8* strType;
    21802180        PS_TYPE_NAME(strType,x->type.type);
    21812181        psError(PS_ERR_BAD_PARAMETER_TYPE,
  • trunk/psLib/src/math/psPolynomial.h

    r2600 r2788  
    1212*  @author GLG, MHPCC
    1313*
    14 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-12-02 21:12:51 $
     14*  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-12-22 05:09:32 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3939 *        \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    4040 *
    41  *  @return float      value on the gaussian curve given the input parameters
    42  */
    43 float psGaussian(
    44     float x,                           ///< Value at which to evaluate
    45     float mean,                        ///< Mean for the Gaussian
    46     float stddev,                      ///< Standard deviation for the Gaussian
     41 *  @return psF32      value on the gaussian curve given the input parameters
     42 */
     43psF32 psGaussian(
     44    psF32 x,                           ///< Value at which to evaluate
     45    psF32 mean,                        ///< Mean for the Gaussian
     46    psF32 stddev,                      ///< Standard deviation for the Gaussian
    4747    psBool normal                        ///< Indicates whether result should be normalized
    4848);
     
    5555 */
    5656psVector* p_psGaussianDev(
    57     float mean,                        ///< The mean of the Gaussian
    58     float sigma,                       ///< The sigma of the Gaussian
     57    psF32 mean,                        ///< The mean of the Gaussian
     58    psF32 sigma,                       ///< The sigma of the Gaussian
    5959    psS32 Npts                           ///< The size of the vector
    6060);
     
    7070    psPolynomialType type;             ///< Polynomial type
    7171    psS32 n;                             ///< Number of terms
    72     float *coeff;                      ///< Coefficients
    73     float *coeffErr;                   ///< Error in coefficients
    74     char *mask;                        ///< Coefficient mask
     72    psF32 *coeff;                      ///< Coefficients
     73    psF32 *coeffErr;                   ///< Error in coefficients
     74    psU8 *mask;                        ///< Coefficient mask
    7575}
    7676psPolynomial1D;
     
    8282    psS32 nX;                            ///< Number of terms in x
    8383    psS32 nY;                            ///< Number of terms in y
    84     float **coeff;                     ///< Coefficients
    85     float **coeffErr;                  ///< Error in coefficients
    86     char **mask;                       ///< Coefficients mask
     84    psF32 **coeff;                     ///< Coefficients
     85    psF32 **coeffErr;                  ///< Error in coefficients
     86    psU8 **mask;                       ///< Coefficients mask
    8787}
    8888psPolynomial2D;
     
    9595    psS32 nY;                            ///< Number of terms in y
    9696    psS32 nZ;                            ///< Number of terms in z
    97     float ***coeff;                    ///< Coefficients
    98     float ***coeffErr;                 ///< Error in coefficients
    99     char ***mask;                      ///< Coefficients mask
     97    psF32 ***coeff;                    ///< Coefficients
     98    psF32 ***coeffErr;                 ///< Error in coefficients
     99    psU8 ***mask;                      ///< Coefficients mask
    100100}
    101101psPolynomial3D;
     
    109109    psS32 nY;                            ///< Number of terms in y
    110110    psS32 nZ;                            ///< Number of terms in z
    111     float ****coeff;                   ///< Coefficients
    112     float ****coeffErr;                ///< Error in coefficients
    113     char ****mask;                     ///< Coefficients mask
     111    psF32 ****coeff;                   ///< Coefficients
     112    psF32 ****coeffErr;                ///< Error in coefficients
     113    psU8 ****mask;                     ///< Coefficients mask
    114114}
    115115psPolynomial4D;
     
    160160/** Evaluates a 1-D polynomial at specific coordinates.
    161161 *
    162  *  @return float    result of polynomial at given location
    163  */
    164 float psPolynomial1DEval(
     162 *  @return psF32    result of polynomial at given location
     163 */
     164psF32 psPolynomial1DEval(
    165165    const psPolynomial1D* myPoly,       ///< Coefficients for the polynomial
    166     float x                           ///< location at which to evaluate
     166    psF32 x                           ///< location at which to evaluate
    167167);
    168168
    169169/** Evaluates a 2-D polynomial at specific coordinates.
    170170 *
    171  *  @return float    result of polynomial at given location
    172  */
    173 float psPolynomial2DEval(
     171 *  @return psF32    result of polynomial at given location
     172 */
     173psF32 psPolynomial2DEval(
    174174    const psPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    175     float x,                           ///< x location at which to evaluate
    176     float y                           ///< y location at which to evaluate
     175    psF32 x,                           ///< x location at which to evaluate
     176    psF32 y                           ///< y location at which to evaluate
    177177);
    178178
    179179/** Evaluates a 3-D polynomial at specific coordinates.
    180180 *
    181  *  @return float    result of polynomial at given location
    182  */
    183 float psPolynomial3DEval(
     181 *  @return psF32    result of polynomial at given location
     182 */
     183psF32 psPolynomial3DEval(
    184184    const psPolynomial3D* myPoly,       ///< Coefficients for the polynomial
    185     float x,                           ///< x location at which to evaluate
    186     float y,                           ///< y location at which to evaluate
    187     float z                           ///< z location at which to evaluate
     185    psF32 x,                           ///< x location at which to evaluate
     186    psF32 y,                           ///< y location at which to evaluate
     187    psF32 z                           ///< z location at which to evaluate
    188188);
    189189
    190190/** Evaluates a 4-D polynomial at specific coordinates.
    191191 *
    192  *  @return float    result of polynomial at given location
    193  */
    194 float psPolynomial4DEval(
     192 *  @return psF32    result of polynomial at given location
     193 */
     194psF32 psPolynomial4DEval(
    195195    const psPolynomial4D* myPoly,       ///< Coefficients for the polynomial
    196     float w,                           ///< w location at which to evaluate
    197     float x,                           ///< x location at which to evaluate
    198     float y,                           ///< y location at which to evaluate
    199     float z                           ///< z location at which to evaluate
     196    psF32 w,                           ///< w location at which to evaluate
     197    psF32 x,                           ///< x location at which to evaluate
     198    psF32 y,                           ///< y location at which to evaluate
     199    psF32 z                           ///< z location at which to evaluate
    200200);
    201201
     
    235235    psPolynomialType type;             ///< Polynomial type
    236236    psS32 n;                             ///< Number of terms
    237     double *coeff;                     ///< Coefficients
    238     double *coeffErr;                  ///< Error in coefficients
    239     char *mask;                        ///< Coefficient mask
     237    psF64 *coeff;                     ///< Coefficients
     238    psF64 *coeffErr;                  ///< Error in coefficients
     239    psU8 *mask;                        ///< Coefficient mask
    240240}
    241241psDPolynomial1D;
     
    247247    psS32 nX;                            ///< Number of terms in x
    248248    psS32 nY;                            ///< Number of terms in y
    249     double **coeff;                    ///< Coefficients
    250     double **coeffErr;                 ///< Error in coefficients
    251     char **mask;                       ///< Coefficients mask
     249    psF64 **coeff;                    ///< Coefficients
     250    psF64 **coeffErr;                 ///< Error in coefficients
     251    psU8 **mask;                       ///< Coefficients mask
    252252}
    253253psDPolynomial2D;
     
    260260    psS32 nY;                            ///< Number of terms in y
    261261    psS32 nZ;                            ///< Number of terms in z
    262     double ***coeff;                   ///< Coefficients
    263     double ***coeffErr;                ///< Error in coefficients
    264     char ***mask;                      ///< Coefficient mask
     262    psF64 ***coeff;                   ///< Coefficients
     263    psF64 ***coeffErr;                ///< Error in coefficients
     264    psU8 ***mask;                      ///< Coefficient mask
    265265}
    266266psDPolynomial3D;
     
    274274    psS32 nY;                            ///< Number of terms in y
    275275    psS32 nZ;                            ///< Number of terms in z
    276     double ****coeff;                  ///< Coefficients
    277     double ****coeffErr;               ///< Error in coefficients
    278     char ****mask;                     ///< Coefficients mask
     276    psF64 ****coeff;                  ///< Coefficients
     277    psF64 ****coeffErr;               ///< Error in coefficients
     278    psU8 ****mask;                     ///< Coefficients mask
    279279}
    280280psDPolynomial4D;
     
    324324/** Evaluates a double-precision 1-D polynomial at specific coordinates.
    325325 *
    326  *  @return float    result of polynomial at given location
    327  */
    328 double psDPolynomial1DEval(
     326 *  @return psF32    result of polynomial at given location
     327 */
     328psF64 psDPolynomial1DEval(
    329329    const psDPolynomial1D* myPoly,      ///< Coefficients for the polynomial
    330     double x                          ///< Value at which to evaluate
     330    psF64 x                          ///< Value at which to evaluate
    331331);
    332332
    333333/** Evaluates a double-precision 2-D polynomial at specific coordinates.
    334334 *
    335  *  @return float    result of polynomial at given location
    336  */
    337 double psDPolynomial2DEval(
     335 *  @return psF32    result of polynomial at given location
     336 */
     337psF64 psDPolynomial2DEval(
    338338    const psDPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    339     double x,                           ///< Value x at which to evaluate
    340     double y            ///< Value y at which to evaluate
     339    psF64 x,                           ///< Value x at which to evaluate
     340    psF64 y            ///< Value y at which to evaluate
    341341);
    342342
    343343/** Evaluates a double-precision 3-D polynomial at specific coordinates.
    344344 *
    345  *  @return float    result of polynomial at given location
    346  */
    347 double psDPolynomial3DEval(
     345 *  @return psF32    result of polynomial at given location
     346 */
     347psF64 psDPolynomial3DEval(
    348348    const psDPolynomial3D* myPoly,      ///< Coefficients for the polynomial
    349     double x,                          ///< Value x at which to evaluate
    350     double y,                          ///< Value y at which to evaluate
    351     double z     ///< Value z at which to evaluate
     349    psF64 x,                          ///< Value x at which to evaluate
     350    psF64 y,                          ///< Value y at which to evaluate
     351    psF64 z     ///< Value z at which to evaluate
    352352);
    353353
    354354/** Evaluates a double-precision 4-D polynomial at specific coordinates.
    355355 *
    356  *  @return float    result of polynomial at given location
    357  */
    358 double psDPolynomial4DEval(
     356 *  @return psF32    result of polynomial at given location
     357 */
     358psF64 psDPolynomial4DEval(
    359359    const psDPolynomial4D* myPoly,      ///< Coefficients for the polynomial
    360     double w,                          ///< Value w at which to evaluate
    361     double x,                          ///< Value x at which to evaluate
    362     double y,                          ///< Value y at which to evaluate
    363     double z     ///< Value z at which to evaluate
     360    psF64 w,                          ///< Value w at which to evaluate
     361    psF64 x,                          ///< Value x at which to evaluate
     362    psF64 y,                          ///< Value y at which to evaluate
     363    psF64 z     ///< Value z at which to evaluate
    364364);
    365365
     
    396396    psS32 n;                        ///< The number of spline polynomials
    397397    psPolynomial1D **spline;      ///< An array of n pointers to the spline polynomials
    398     float *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    399     float *domains;               ///< The boundaries between each spline piece.  Size is n+1.
     398    psF32 *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
     399    psF32 *domains;               ///< The boundaries between each spline piece.  Size is n+1.
    400400}
    401401psSpline1D;
     
    403403psSpline1D *psSpline1DAlloc(psS32 n,
    404404                            psS32 order,
    405                             float min,
    406                             float max);
     405                            psF32 min,
     406                            psF32 max);
    407407
    408408psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
    409409                                   psS32 order);
    410410
    411 float psSpline1DEval(
    412     float x,
     411psF32 psSpline1DEval(
     412    psF32 x,
    413413    const psSpline1D *spline
    414414);
     
    427427                                psScalar *x);
    428428
    429 float p_psNRSpline1DEval(psSpline1D *spline,
     429psF32 p_psNRSpline1DEval(psSpline1D *spline,
    430430                         const psVector* restrict x,
    431431                         const psVector* restrict y,
    432                          float X);
     432                         psF32 X);
    433433
    434434/* \} */// End of MathGroup Functions
  • trunk/psLib/src/math/psSpline.c

    r2778 r2788  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-12-21 20:42:07 $
     9 *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-12-22 05:09:32 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5555static void dPolynomial4DFree(psDPolynomial4D* myPoly);
    5656static void spline1DFree(psSpline1D *tmpSpline);
    57 static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x);
     57static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
    5858static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
    5959
     
    277277 
    278278 *****************************************************************************/
    279 static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     279static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    280280{
    281281    psS32 loop_x = 0;
    282     float polySum = 0.0;
    283     float xSum = 1.0;
     282    psF32 polySum = 0.0;
     283    psF32 xSum = 1.0;
    284284
    285285    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
     
    307307// XXX: How does the mask vector effect Crenshaw's formula?
    308308// XXX: We assume that x is scaled between -1.0 and 1.0;
    309 static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
     309static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
    310310{
    311311    //    PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0);
     
    314314    psS32 n;
    315315    psS32 i;
    316     float tmp;
     316    psF32 tmp;
    317317
    318318    n = myPoly->n;
     
    337337    psS32 n;
    338338    psS32 i;
    339     float tmp;
     339    psF32 tmp;
    340340    psPolynomial1D **chebPolys = NULL;
    341341
     
    355355}
    356356
    357 static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     357static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    358358{
    359359    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    361361    psS32 loop_x = 0;
    362362    psS32 loop_y = 0;
    363     float polySum = 0.0;
    364     float xSum = 1.0;
    365     float ySum = 1.0;
     363    psF32 polySum = 0.0;
     364    psF32 xSum = 1.0;
     365    psF32 ySum = 1.0;
    366366
    367367    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    379379}
    380380
    381 static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
     381static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
    382382{
    383383    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    386386    psS32 loop_y = 0;
    387387    psS32 i = 0;
    388     float polySum = 0.0;
     388    psF32 polySum = 0.0;
    389389    psPolynomial1D* *chebPolys = NULL;
    390390    psS32 maxChebyPoly = 0;
     
    414414}
    415415
    416 static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     416static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    417417{
    418418    psS32 loop_x = 0;
    419419    psS32 loop_y = 0;
    420420    psS32 loop_z = 0;
    421     float polySum = 0.0;
    422     float xSum = 1.0;
    423     float ySum = 1.0;
    424     float zSum = 1.0;
     421    psF32 polySum = 0.0;
     422    psF32 xSum = 1.0;
     423    psF32 ySum = 1.0;
     424    psF32 zSum = 1.0;
    425425
    426426    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    442442}
    443443
    444 static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
     444static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
    445445{
    446446    psS32 loop_x = 0;
     
    448448    psS32 loop_z = 0;
    449449    psS32 i = 0;
    450     float polySum = 0.0;
     450    psF32 polySum = 0.0;
    451451    psPolynomial1D* *chebPolys = NULL;
    452452    psS32 maxChebyPoly = 0;
     
    483483}
    484484
    485 static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     485static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    486486{
    487487    psS32 loop_w = 0;
     
    489489    psS32 loop_y = 0;
    490490    psS32 loop_z = 0;
    491     float polySum = 0.0;
    492     float wSum = 1.0;
    493     float xSum = 1.0;
    494     float ySum = 1.0;
    495     float zSum = 1.0;
     491    psF32 polySum = 0.0;
     492    psF32 wSum = 1.0;
     493    psF32 xSum = 1.0;
     494    psF32 ySum = 1.0;
     495    psF32 zSum = 1.0;
    496496
    497497    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    517517}
    518518
    519 static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
     519static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
    520520{
    521521    psS32 loop_w = 0;
     
    524524    psS32 loop_z = 0;
    525525    psS32 i = 0;
    526     float polySum = 0.0;
     526    psF32 polySum = 0.0;
    527527    psPolynomial1D* *chebPolys = NULL;
    528528    psS32 maxChebyPoly = 0;
     
    568568    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    569569 *****************************************************************************/
    570 static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     570static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    571571{
    572572    psS32 loop_x = 0;
    573     double polySum = 0.0;
    574     double xSum = 1.0;
     573    psF64 polySum = 0.0;
     574    psF64 xSum = 1.0;
    575575
    576576    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     
    586586// XXX: You can do this without having to psAlloc() vector d.
    587587// XXX: How does the mask vector effect Crenshaw's formula?
    588 static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
     588static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
    589589{
    590590    psVector *d;
    591591    psS32 n;
    592592    psS32 i;
    593     double tmp;
     593    psF64 tmp;
    594594
    595595    n = myPoly->n;
     
    611611}
    612612
    613 static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     613static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    614614{
    615615    psS32 loop_x = 0;
    616616    psS32 loop_y = 0;
    617     double polySum = 0.0;
    618     double xSum = 1.0;
    619     double ySum = 1.0;
     617    psF64 polySum = 0.0;
     618    psF64 xSum = 1.0;
     619    psF64 ySum = 1.0;
    620620
    621621    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    633633}
    634634
    635 static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
     635static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
    636636{
    637637    psS32 loop_x = 0;
    638638    psS32 loop_y = 0;
    639639    psS32 i = 0;
    640     double polySum = 0.0;
     640    psF64 polySum = 0.0;
    641641    psPolynomial1D* *chebPolys = NULL;
    642642    psS32 maxChebyPoly = 0;
     
    667667}
    668668
    669 static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     669static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    670670{
    671671    psS32 loop_x = 0;
    672672    psS32 loop_y = 0;
    673673    psS32 loop_z = 0;
    674     double polySum = 0.0;
    675     double xSum = 1.0;
    676     double ySum = 1.0;
    677     double zSum = 1.0;
     674    psF64 polySum = 0.0;
     675    psF64 xSum = 1.0;
     676    psF64 ySum = 1.0;
     677    psF64 zSum = 1.0;
    678678
    679679    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
     
    695695}
    696696
    697 static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
     697static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
    698698{
    699699    psS32 loop_x = 0;
     
    701701    psS32 loop_z = 0;
    702702    psS32 i = 0;
    703     double polySum = 0.0;
     703    psF64 polySum = 0.0;
    704704    psPolynomial1D* *chebPolys = NULL;
    705705    psS32 maxChebyPoly = 0;
     
    736736}
    737737
    738 static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     738static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    739739{
    740740    psS32 loop_w = 0;
     
    742742    psS32 loop_y = 0;
    743743    psS32 loop_z = 0;
    744     double polySum = 0.0;
    745     double wSum = 1.0;
    746     double xSum = 1.0;
    747     double ySum = 1.0;
    748     double zSum = 1.0;
     744    psF64 polySum = 0.0;
     745    psF64 wSum = 1.0;
     746    psF64 xSum = 1.0;
     747    psF64 ySum = 1.0;
     748    psF64 zSum = 1.0;
    749749
    750750    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
     
    770770}
    771771
    772 static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
     772static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
    773773{
    774774    psS32 loop_w = 0;
     
    777777    psS32 loop_z = 0;
    778778    psS32 i = 0;
    779     double polySum = 0.0;
     779    psF64 polySum = 0.0;
    780780    psPolynomial1D* *chebPolys = NULL;
    781781    psS32 maxChebyPoly = 0;
     
    828828 *****************************************************************************/
    829829#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
    830 static float fullInterpolate1D##TYPE(ps##TYPE *domain, \
     830static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
    831831                                     ps##TYPE *range, \
    832832                                     psS32 n, \
     
    896896LaGrange interpolation.
    897897 *****************************************************************************/
    898 static float interpolate1DF32(float *domain,
    899                               float *range,
     898static psF32 interpolate1DF32(psF32 *domain,
     899                              psF32 *range,
    900900                              psS32 n,
    901901                              psS32 order,
    902                               float x)
     902                              psF32 x)
    903903{
    904904    psS32 binNum;
     
    941941    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    942942 *****************************************************************************/
    943 float psGaussian(float x, float mean, float sigma, psBool normal)
    944 {
    945     float tmp = 1.0;
     943psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal)
     944{
     945    psF32 tmp = 1.0;
    946946
    947947    psTrace(".psLib.dataManip.psFunctions.psGaussian", 4,
     
    967967XXX: There is no way to seed the random generator.
    968968 *****************************************************************************/
    969 psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts)
     969psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
    970970{
    971971    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
     
    10071007    newPoly->type = type;
    10081008    newPoly->n = n;
    1009     newPoly->coeff = (float *)psAlloc(n * sizeof(float));
    1010     newPoly->coeffErr = (float *)psAlloc(n * sizeof(float));
    1011     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1009    newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32));
     1010    newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32));
     1011    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    10121012    for (i = 0; i < n; i++) {
    10131013        newPoly->coeff[i] = 0.0;
     
    10361036    newPoly->nY = nY;
    10371037
    1038     newPoly->coeff = (float **)psAlloc(nX * sizeof(float *));
    1039     newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *));
    1040     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1038    newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1039    newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     1040    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    10411041    for (x = 0; x < nX; x++) {
    1042         newPoly->coeff[x] = (float *)psAlloc(nY * sizeof(float));
    1043         newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float));
    1044         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1042        newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1043        newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     1044        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    10451045    }
    10461046    for (x = 0; x < nX; x++) {
     
    10751075    newPoly->nZ = nZ;
    10761076
    1077     newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **));
    1078     newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **));
    1079     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1077    newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1078    newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1079    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    10801080    for (x = 0; x < nX; x++) {
    1081         newPoly->coeff[x] = (float **)psAlloc(nY * sizeof(float *));
    1082         newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *));
    1083         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1081        newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1082        newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1083        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    10841084        for (y = 0; y < nY; y++) {
    1085             newPoly->coeff[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1086             newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float));
    1087             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1085            newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1086            newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1087            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    10881088        }
    10891089    }
     
    11241124    newPoly->nZ = nZ;
    11251125
    1126     newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***));
    1127     newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***));
    1128     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1126    newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1127    newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
     1128    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    11291129    for (w = 0; w < nW; w++) {
    1130         newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **));
    1131         newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **));
    1132         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1130        newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1131        newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     1132        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    11331133        for (x = 0; x < nX; x++) {
    1134             newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1135             newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *));
    1136             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1134            newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1135            newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     1136            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    11371137            for (y = 0; y < nY; y++) {
    1138                 newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1139                 newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
    1140                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1138                newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1139                newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     1140                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    11411141            }
    11421142        }
     
    11571157}
    11581158
    1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, float x)
     1159psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x)
    11601160{
    11611161    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    11951195}
    11961196
    1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y)
     1197psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y)
    11981198{
    11991199    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    12501250}
    12511251
    1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z)
     1252psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z)
    12531253{
    12541254    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    13161316}
    13171317
    1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z)
     1318psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z)
    13191319{
    13201320    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    14061406    newPoly->type = type;
    14071407    newPoly->n = n;
    1408     newPoly->coeff = (double *)psAlloc(n * sizeof(double));
    1409     newPoly->coeffErr = (double *)psAlloc(n * sizeof(double));
    1410     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     1408    newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64));
     1409    newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64));
     1410    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
    14111411    for (i = 0; i < n; i++) {
    14121412        newPoly->coeff[i] = 0.0;
     
    14351435    newPoly->nY = nY;
    14361436
    1437     newPoly->coeff = (double **)psAlloc(nX * sizeof(double *));
    1438     newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *));
    1439     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     1437    newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1438    newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *));
     1439    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
    14401440    for (x = 0; x < nX; x++) {
    1441         newPoly->coeff[x] = (double *)psAlloc(nY * sizeof(double));
    1442         newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double));
    1443         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     1441        newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1442        newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
     1443        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
    14441444    }
    14451445    for (x = 0; x < nX; x++) {
     
    14741474    newPoly->nZ = nZ;
    14751475
    1476     newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **));
    1477     newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **));
    1478     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     1476    newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1477    newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1478    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    14791479    for (x = 0; x < nX; x++) {
    1480         newPoly->coeff[x] = (double **)psAlloc(nY * sizeof(double *));
    1481         newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *));
    1482         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     1480        newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1481        newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1482        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    14831483        for (y = 0; y < nY; y++) {
    1484             newPoly->coeff[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1485             newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double));
    1486             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     1484            newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1485            newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1486            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    14871487        }
    14881488    }
     
    15231523    newPoly->nZ = nZ;
    15241524
    1525     newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***));
    1526     newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***));
    1527     newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
     1525    newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1526    newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
     1527    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
    15281528    for (w = 0; w < nW; w++) {
    1529         newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **));
    1530         newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **));
    1531         newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
     1529        newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1530        newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
     1531        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
    15321532        for (x = 0; x < nX; x++) {
    1533             newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1534             newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *));
    1535             newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
     1533            newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1534            newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
     1535            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
    15361536            for (y = 0; y < nY; y++) {
    1537                 newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1538                 newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
    1539                 newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
     1537                newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1538                newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
     1539                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
    15401540            }
    15411541        }
     
    15571557
    15581558
    1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x)
     1559psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x)
    15601560{
    15611561    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    15981598
    15991599
    1600 double psDPolynomial2DEval(const psDPolynomial2D* myPoly,
    1601                            double x,
    1602                            double y)
     1600psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly,
     1601                          psF64 x,
     1602                          psF64 y)
    16031603{
    16041604    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    16551655
    16561656
    1657 double psDPolynomial3DEval(const psDPolynomial3D* myPoly,
    1658                            double x,
    1659                            double y,
    1660                            double z)
     1657psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly,
     1658                          psF64 x,
     1659                          psF64 y,
     1660                          psF64 z)
    16611661{
    16621662    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    17241724}
    17251725
    1726 double psDPolynomial4DEval(const psDPolynomial4D* myPoly,
    1727                            double w,
    1728                            double x,
    1729                            double y,
    1730                            double z)
     1726psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly,
     1727                          psF64 w,
     1728                          psF64 x,
     1729                          psF64 y,
     1730                          psF64 z)
    17311731{
    17321732    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    18101810//    psS32 n;
    18111811//    psPolynomial1D **spline;
    1812 //    float *p_psDeriv2;
    1813 //    float *domains;
     1812//    psF32 *p_psDeriv2;
     1813//    psF32 *domains;
    18141814//} psSpline1D;
    18151815
     
    18221822psSpline1D *psSpline1DAlloc(psS32 numSplines,
    18231823                            psS32 order,
    1824                             float min,
    1825                             float max)
     1824                            psF32 min,
     1825                            psF32 max)
    18261826{
    18271827    PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL);
     
    18311831    psSpline1D *tmp = NULL;
    18321832    psS32 i;
    1833     float tmpDomain;
    1834     float width;
     1833    psF32 tmpDomain;
     1834    psF32 width;
    18351835
    18361836    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
     
    18451845    tmp->p_psDeriv2 = NULL;
    18461846
    1847     tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    1848     width = (max - min) / ((float) numSplines);
     1847    tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32));
     1848    width = (max - min) / ((psF32) numSplines);
    18491849
    18501850    (tmp->domains)[0] = min;
     
    18841884    }
    18851885
    1886     tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float));
     1886    tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));
    18871887
    18881888    for (i=0;i<bounds->n;i++) {
     
    19771977{
    19781978    PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3);
    1979     char* strType;
     1979    psS8* strType;
    19801980
    19811981    switch (x->type.type) {
     
    20742074        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
    20752075
    2076         psScalar *tmpScalar = psScalarAlloc((double)
     2076        psScalar *tmpScalar = psScalarAlloc((psF64)
    20772077                                            interpolate1DF32(domain32->data.F32,
    20782078                                                             range32->data.F32,
    20792079                                                             domain32->n,
    20802080                                                             order,
    2081                                                              (float) x->data.F64), PS_TYPE_F64);
     2081                                                             (psF32) x->data.F64), PS_TYPE_F64);
    20822082        psFree(range32);
    20832083        psFree(domain32);
     
    20892089
    20902090    } else {
    2091         char* strType;
     2091        psS8* strType;
    20922092        PS_TYPE_NAME(strType,x->type.type);
    20932093        psError(PS_ERR_BAD_PARAMETER_TYPE,
     
    21152115     the spline fit functions require F32 and F64.
    21162116 *****************************************************************************/
    2117 float psSpline1DEval(
    2118     float x,
     2117psF32 psSpline1DEval(
     2118    psF32 x,
    21192119    const psSpline1D *spline
    21202120)
     
    21722172        for (i=0;i<x->n;i++) {
    21732173            tmpVector->data.F32[i] = psSpline1DEval(
    2174                                          (float) x->data.F64[i],
     2174                                         (psF32) x->data.F64[i],
    21752175                                         spline
    21762176                                     );
    21772177        }
    21782178    } else {
    2179         char* strType;
     2179        psS8* strType;
    21802180        PS_TYPE_NAME(strType,x->type.type);
    21812181        psError(PS_ERR_BAD_PARAMETER_TYPE,
  • trunk/psLib/src/math/psSpline.h

    r2600 r2788  
    1212*  @author GLG, MHPCC
    1313*
    14 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-12-02 21:12:51 $
     14*  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-12-22 05:09:32 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3939 *        \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    4040 *
    41  *  @return float      value on the gaussian curve given the input parameters
    42  */
    43 float psGaussian(
    44     float x,                           ///< Value at which to evaluate
    45     float mean,                        ///< Mean for the Gaussian
    46     float stddev,                      ///< Standard deviation for the Gaussian
     41 *  @return psF32      value on the gaussian curve given the input parameters
     42 */
     43psF32 psGaussian(
     44    psF32 x,                           ///< Value at which to evaluate
     45    psF32 mean,                        ///< Mean for the Gaussian
     46    psF32 stddev,                      ///< Standard deviation for the Gaussian
    4747    psBool normal                        ///< Indicates whether result should be normalized
    4848);
     
    5555 */
    5656psVector* p_psGaussianDev(
    57     float mean,                        ///< The mean of the Gaussian
    58     float sigma,                       ///< The sigma of the Gaussian
     57    psF32 mean,                        ///< The mean of the Gaussian
     58    psF32 sigma,                       ///< The sigma of the Gaussian
    5959    psS32 Npts                           ///< The size of the vector
    6060);
     
    7070    psPolynomialType type;             ///< Polynomial type
    7171    psS32 n;                             ///< Number of terms
    72     float *coeff;                      ///< Coefficients
    73     float *coeffErr;                   ///< Error in coefficients
    74     char *mask;                        ///< Coefficient mask
     72    psF32 *coeff;                      ///< Coefficients
     73    psF32 *coeffErr;                   ///< Error in coefficients
     74    psU8 *mask;                        ///< Coefficient mask
    7575}
    7676psPolynomial1D;
     
    8282    psS32 nX;                            ///< Number of terms in x
    8383    psS32 nY;                            ///< Number of terms in y
    84     float **coeff;                     ///< Coefficients
    85     float **coeffErr;                  ///< Error in coefficients
    86     char **mask;                       ///< Coefficients mask
     84    psF32 **coeff;                     ///< Coefficients
     85    psF32 **coeffErr;                  ///< Error in coefficients
     86    psU8 **mask;                       ///< Coefficients mask
    8787}
    8888psPolynomial2D;
     
    9595    psS32 nY;                            ///< Number of terms in y
    9696    psS32 nZ;                            ///< Number of terms in z
    97     float ***coeff;                    ///< Coefficients
    98     float ***coeffErr;                 ///< Error in coefficients
    99     char ***mask;                      ///< Coefficients mask
     97    psF32 ***coeff;                    ///< Coefficients
     98    psF32 ***coeffErr;                 ///< Error in coefficients
     99    psU8 ***mask;                      ///< Coefficients mask
    100100}
    101101psPolynomial3D;
     
    109109    psS32 nY;                            ///< Number of terms in y
    110110    psS32 nZ;                            ///< Number of terms in z
    111     float ****coeff;                   ///< Coefficients
    112     float ****coeffErr;                ///< Error in coefficients
    113     char ****mask;                     ///< Coefficients mask
     111    psF32 ****coeff;                   ///< Coefficients
     112    psF32 ****coeffErr;                ///< Error in coefficients
     113    psU8 ****mask;                     ///< Coefficients mask
    114114}
    115115psPolynomial4D;
     
    160160/** Evaluates a 1-D polynomial at specific coordinates.
    161161 *
    162  *  @return float    result of polynomial at given location
    163  */
    164 float psPolynomial1DEval(
     162 *  @return psF32    result of polynomial at given location
     163 */
     164psF32 psPolynomial1DEval(
    165165    const psPolynomial1D* myPoly,       ///< Coefficients for the polynomial
    166     float x                           ///< location at which to evaluate
     166    psF32 x                           ///< location at which to evaluate
    167167);
    168168
    169169/** Evaluates a 2-D polynomial at specific coordinates.
    170170 *
    171  *  @return float    result of polynomial at given location
    172  */
    173 float psPolynomial2DEval(
     171 *  @return psF32    result of polynomial at given location
     172 */
     173psF32 psPolynomial2DEval(
    174174    const psPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    175     float x,                           ///< x location at which to evaluate
    176     float y                           ///< y location at which to evaluate
     175    psF32 x,                           ///< x location at which to evaluate
     176    psF32 y                           ///< y location at which to evaluate
    177177);
    178178
    179179/** Evaluates a 3-D polynomial at specific coordinates.
    180180 *
    181  *  @return float    result of polynomial at given location
    182  */
    183 float psPolynomial3DEval(
     181 *  @return psF32    result of polynomial at given location
     182 */
     183psF32 psPolynomial3DEval(
    184184    const psPolynomial3D* myPoly,       ///< Coefficients for the polynomial
    185     float x,                           ///< x location at which to evaluate
    186     float y,                           ///< y location at which to evaluate
    187     float z                           ///< z location at which to evaluate
     185    psF32 x,                           ///< x location at which to evaluate
     186    psF32 y,                           ///< y location at which to evaluate
     187    psF32 z                           ///< z location at which to evaluate
    188188);
    189189
    190190/** Evaluates a 4-D polynomial at specific coordinates.
    191191 *
    192  *  @return float    result of polynomial at given location
    193  */
    194 float psPolynomial4DEval(
     192 *  @return psF32    result of polynomial at given location
     193 */
     194psF32 psPolynomial4DEval(
    195195    const psPolynomial4D* myPoly,       ///< Coefficients for the polynomial
    196     float w,                           ///< w location at which to evaluate
    197     float x,                           ///< x location at which to evaluate
    198     float y,                           ///< y location at which to evaluate
    199     float z                           ///< z location at which to evaluate
     196    psF32 w,                           ///< w location at which to evaluate
     197    psF32 x,                           ///< x location at which to evaluate
     198    psF32 y,                           ///< y location at which to evaluate
     199    psF32 z                           ///< z location at which to evaluate
    200200);
    201201
     
    235235    psPolynomialType type;             ///< Polynomial type
    236236    psS32 n;                             ///< Number of terms
    237     double *coeff;                     ///< Coefficients
    238     double *coeffErr;                  ///< Error in coefficients
    239     char *mask;                        ///< Coefficient mask
     237    psF64 *coeff;                     ///< Coefficients
     238    psF64 *coeffErr;                  ///< Error in coefficients
     239    psU8 *mask;                        ///< Coefficient mask
    240240}
    241241psDPolynomial1D;
     
    247247    psS32 nX;                            ///< Number of terms in x
    248248    psS32 nY;                            ///< Number of terms in y
    249     double **coeff;                    ///< Coefficients
    250     double **coeffErr;                 ///< Error in coefficients
    251     char **mask;                       ///< Coefficients mask
     249    psF64 **coeff;                    ///< Coefficients
     250    psF64 **coeffErr;                 ///< Error in coefficients
     251    psU8 **mask;                       ///< Coefficients mask
    252252}
    253253psDPolynomial2D;
     
    260260    psS32 nY;                            ///< Number of terms in y
    261261    psS32 nZ;                            ///< Number of terms in z
    262     double ***coeff;                   ///< Coefficients
    263     double ***coeffErr;                ///< Error in coefficients
    264     char ***mask;                      ///< Coefficient mask
     262    psF64 ***coeff;                   ///< Coefficients
     263    psF64 ***coeffErr;                ///< Error in coefficients
     264    psU8 ***mask;                      ///< Coefficient mask
    265265}
    266266psDPolynomial3D;
     
    274274    psS32 nY;                            ///< Number of terms in y
    275275    psS32 nZ;                            ///< Number of terms in z
    276     double ****coeff;                  ///< Coefficients
    277     double ****coeffErr;               ///< Error in coefficients
    278     char ****mask;                     ///< Coefficients mask
     276    psF64 ****coeff;                  ///< Coefficients
     277    psF64 ****coeffErr;               ///< Error in coefficients
     278    psU8 ****mask;                     ///< Coefficients mask
    279279}
    280280psDPolynomial4D;
     
    324324/** Evaluates a double-precision 1-D polynomial at specific coordinates.
    325325 *
    326  *  @return float    result of polynomial at given location
    327  */
    328 double psDPolynomial1DEval(
     326 *  @return psF32    result of polynomial at given location
     327 */
     328psF64 psDPolynomial1DEval(
    329329    const psDPolynomial1D* myPoly,      ///< Coefficients for the polynomial
    330     double x                          ///< Value at which to evaluate
     330    psF64 x                          ///< Value at which to evaluate
    331331);
    332332
    333333/** Evaluates a double-precision 2-D polynomial at specific coordinates.
    334334 *
    335  *  @return float    result of polynomial at given location
    336  */
    337 double psDPolynomial2DEval(
     335 *  @return psF32    result of polynomial at given location
     336 */
     337psF64 psDPolynomial2DEval(
    338338    const psDPolynomial2D* myPoly,       ///< Coefficients for the polynomial
    339     double x,                           ///< Value x at which to evaluate
    340     double y            ///< Value y at which to evaluate
     339    psF64 x,                           ///< Value x at which to evaluate
     340    psF64 y            ///< Value y at which to evaluate
    341341);
    342342
    343343/** Evaluates a double-precision 3-D polynomial at specific coordinates.
    344344 *
    345  *  @return float    result of polynomial at given location
    346  */
    347 double psDPolynomial3DEval(
     345 *  @return psF32    result of polynomial at given location
     346 */
     347psF64 psDPolynomial3DEval(
    348348    const psDPolynomial3D* myPoly,      ///< Coefficients for the polynomial
    349     double x,                          ///< Value x at which to evaluate
    350     double y,                          ///< Value y at which to evaluate
    351     double z     ///< Value z at which to evaluate
     349    psF64 x,                          ///< Value x at which to evaluate
     350    psF64 y,                          ///< Value y at which to evaluate
     351    psF64 z     ///< Value z at which to evaluate
    352352);
    353353
    354354/** Evaluates a double-precision 4-D polynomial at specific coordinates.
    355355 *
    356  *  @return float    result of polynomial at given location
    357  */
    358 double psDPolynomial4DEval(
     356 *  @return psF32    result of polynomial at given location
     357 */
     358psF64 psDPolynomial4DEval(
    359359    const psDPolynomial4D* myPoly,      ///< Coefficients for the polynomial
    360     double w,                          ///< Value w at which to evaluate
    361     double x,                          ///< Value x at which to evaluate
    362     double y,                          ///< Value y at which to evaluate
    363     double z     ///< Value z at which to evaluate
     360    psF64 w,                          ///< Value w at which to evaluate
     361    psF64 x,                          ///< Value x at which to evaluate
     362    psF64 y,                          ///< Value y at which to evaluate
     363    psF64 z     ///< Value z at which to evaluate
    364364);
    365365
     
    396396    psS32 n;                        ///< The number of spline polynomials
    397397    psPolynomial1D **spline;      ///< An array of n pointers to the spline polynomials
    398     float *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    399     float *domains;               ///< The boundaries between each spline piece.  Size is n+1.
     398    psF32 *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
     399    psF32 *domains;               ///< The boundaries between each spline piece.  Size is n+1.
    400400}
    401401psSpline1D;
     
    403403psSpline1D *psSpline1DAlloc(psS32 n,
    404404                            psS32 order,
    405                             float min,
    406                             float max);
     405                            psF32 min,
     406                            psF32 max);
    407407
    408408psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
    409409                                   psS32 order);
    410410
    411 float psSpline1DEval(
    412     float x,
     411psF32 psSpline1DEval(
     412    psF32 x,
    413413    const psSpline1D *spline
    414414);
     
    427427                                psScalar *x);
    428428
    429 float p_psNRSpline1DEval(psSpline1D *spline,
     429psF32 p_psNRSpline1DEval(psSpline1D *spline,
    430430                         const psVector* restrict x,
    431431                         const psVector* restrict y,
    432                          float X);
     432                         psF32 X);
    433433
    434434/* \} */// End of MathGroup Functions
  • trunk/psLib/src/math/psStats.c

    r2782 r2788  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.107 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-22 00:54:28 $
     11 *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-22 05:09:32 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    7070/*****************************************************************************/
    7171
    72 psBool p_psGetStatValue(const psStats* stats, double *value)
     72psBool p_psGetStatValue(const psStats* stats, psF64 *value)
    7373{
    7474
     
    145145this routine sets stats->sampleMean to NAN.
    146146 *****************************************************************************/
    147 int p_psVectorSampleMean(const psVector* restrict myVector,
    148                          const psVector* restrict errors,
    149                          const psVector* restrict maskVector,
    150                          psU32 maskVal,
    151                          psStats* stats)
     147psS32 p_psVectorSampleMean(const psVector* restrict myVector,
     148                           const psVector* restrict errors,
     149                           const psVector* restrict maskVector,
     150                           psU32 maskVal,
     151                           psStats* stats)
    152152{
    153153
    154154    psS32 i = 0;                // Loop index variable
    155     float mean = 0.0;           // The mean
     155    psF32 mean = 0.0;           // The mean
    156156    psS32 count = 0;            // # of points in this mean
    157157
     
    171171                }
    172172                if (count != 0) {
    173                     mean /= (float)count;
     173                    mean /= (psF32)count;
    174174                } else {
    175175                    mean = NAN;
     
    185185                }
    186186                if (count != 0) {
    187                     mean /= (float)count;
     187                    mean /= (psF32)count;
    188188                } else {
    189189                    mean = NAN;
     
    199199                }
    200200                if (count != 0) {
    201                     mean /= (float)count;
     201                    mean /= (psF32)count;
    202202                } else {
    203203                    mean = NAN;
     
    207207                    mean += myVector->data.F32[i];
    208208                }
    209                 mean /= (float)myVector->n;
     209                mean /= (psF32)myVector->n;
    210210            }
    211211        }
     
    283283this routine sets stats->max to NAN.
    284284 *****************************************************************************/
    285 int p_psVectorMax(const psVector* restrict myVector,
    286                   const psVector* restrict maskVector,
    287                   psU32 maskVal,
    288                   psStats* stats)
     285psS32 p_psVectorMax(const psVector* restrict myVector,
     286                    const psVector* restrict maskVector,
     287                    psU32 maskVal,
     288                    psStats* stats)
    289289{
    290290    psS32 i = 0;                // Loop index variable
    291     float max = -PS_MAX_F32;    // The calculated maximum
     291    psF32 max = -PS_MAX_F32;    // The calculated maximum
    292292    psS32 empty = true;         // Does this vector have valid elements?
    293293
     
    348348this routine sets stats->min to NAN.
    349349 *****************************************************************************/
    350 int p_psVectorMin(const psVector* restrict myVector,
    351                   const psVector* restrict maskVector,
    352                   psU32 maskVal,
    353                   psStats* stats)
     350psS32 p_psVectorMin(const psVector* restrict myVector,
     351                    const psVector* restrict maskVector,
     352                    psU32 maskVal,
     353                    psStats* stats)
    354354{
    355355    psS32 i = 0;                // Loop index variable
    356     float min = PS_MAX_F32;   // The calculated maximum
     356    psF32 min = PS_MAX_F32;   // The calculated maximum
    357357    psS32 empty = true;         // Does this vector have valid elements?
    358358
     
    612612 *****************************************************************************/
    613613psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
    614                                        float sigma)
     614                                       psF32 sigma)
    615615{
    616616    PS_PTR_CHECK_NULL(robustHistogram, NULL);
     
    619619    psS32 i = 0;                  // Loop index variable
    620620    psS32 j = 0;                  // Loop index variable
    621     float iMid;
    622     float jMid;
     621    psF32 iMid;
     622    psF32 jMid;
    623623    psS32 numBins = robustHistogram->nums->n;
    624624    psS32 numBounds = robustHistogram->bounds->n;
     
    626626    psS32 jMin = 0;
    627627    psS32 jMax = 0;
    628     float firstBound = robustHistogram->bounds->data.F32[0];
    629     float lastBound = robustHistogram->bounds->data.F32[numBounds-1];
     628    psF32 firstBound = robustHistogram->bounds->data.F32[0];
     629    psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1];
    630630    psScalar x;
    631631
     
    795795    psS32 i = 0;                  // Loop index variable
    796796    psS32 countInt = 0;           // # of data points being used
    797     float countFloat = 0.0;     // # of data points being used
    798     float mean = 0.0;           // The mean
    799     float diff = 0.0;           // Used in calculating stdev
    800     float sumSquares = 0.0;     // temporary variable
    801     float sumDiffs = 0.0;       // temporary variable
     797    psF32 countFloat = 0.0;     // # of data points being used
     798    psF32 mean = 0.0;           // The mean
     799    psF32 diff = 0.0;           // Used in calculating stdev
     800    psF32 sumSquares = 0.0;     // temporary variable
     801    psF32 sumDiffs = 0.0;       // temporary variable
    802802
    803803    // This procedure requires the mean.  If it has not been already
     
    865865        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
    866866    } else {
    867         countFloat = (float)countInt;
     867        countFloat = (psF32)countInt;
    868868        stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    869869    }
     
    889889    psS32 i = 0;                  // Loop index variable
    890890    psS32 countInt = 0;           // # of data points being used
    891     float countFloat = 0.0;     // # of data points being used
    892     float mean = 0.0;           // The mean
    893     float diff = 0.0;           // Used in calculating stdev
    894     float sumSquares = 0.0;     // temporary variable
    895     float sumDiffs = 0.0;       // temporary variable
     891    psF32 countFloat = 0.0;     // # of data points being used
     892    psF32 mean = 0.0;           // The mean
     893    psF32 diff = 0.0;           // Used in calculating stdev
     894    psF32 sumSquares = 0.0;     // temporary variable
     895    psF32 sumDiffs = 0.0;       // temporary variable
    896896    //    psF32 sum1;
    897897    //    psF32 sum2;
     
    982982            stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor));
    983983        } else {
    984             countFloat = (float)countInt;
     984            countFloat = (psF32)countInt;
    985985            stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    986986
     
    10031003    -2: warning
    10041004 *****************************************************************************/
    1005 int p_psVectorClippedStats(const psVector* restrict myVector,
    1006                            const psVector* restrict errors,
    1007                            const psVector* restrict maskVector,
    1008                            psU32 maskVal,
    1009                            psStats* stats)
     1005psS32 p_psVectorClippedStats(const psVector* restrict myVector,
     1006                             const psVector* restrict errors,
     1007                             const psVector* restrict maskVector,
     1008                             psU32 maskVal,
     1009                             psStats* stats)
    10101010{
    10111011    psS32 i = 0;                  // Loop index variable
    10121012    psS32 j = 0;                  // Loop index variable
    1013     float clippedMean = 0.0;    // self-explanatory
    1014     float clippedStdev = 0.0;   // self-explanatory
    1015     float oldStanMean = 0.0;    // Temporary variable
    1016     float oldStanStdev = 0.0;   // Temporary variable
     1013    psF32 clippedMean = 0.0;    // self-explanatory
     1014    psF32 clippedStdev = 0.0;   // self-explanatory
     1015    psF32 oldStanMean = 0.0;    // Temporary variable
     1016    psF32 oldStanStdev = 0.0;   // Temporary variable
    10171017    psVector* tmpMask = NULL;   // Temporary vector
    10181018
     
    12401240XXX: Create a 2nd-order polynomial version and solve for X analytically.
    12411241 *****************************************************************************/
    1242 float p_ps1DPolyMedian(psPolynomial1D* myPoly,
    1243                        float rangeLow,
    1244                        float rangeHigh,
    1245                        float getThisValue)
     1242psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
     1243                       psF32 rangeLow,
     1244                       psF32 rangeHigh,
     1245                       psF32 getThisValue)
    12461246{
    12471247    PS_POLY_CHECK_NULL(myPoly, NAN);
     
    12501250    // falls within the range of y-values of the polynomial "myPoly" in the
    12511251    // specified x-range (rangeLow:rangeHigh).
    1252     float fLo = psPolynomial1DEval(
     1252    psF32 fLo = psPolynomial1DEval(
    12531253                    myPoly,
    12541254                    rangeLow
    12551255                );
    1256     float fHi = psPolynomial1DEval(
     1256    psF32 fHi = psPolynomial1DEval(
    12571257                    myPoly,
    12581258                    rangeHigh
     
    12661266
    12671267    psS32 numIterations = 0;
    1268     float midpoint = 0.0;
    1269     float oldMidpoint = 1.0;
    1270     float f = 0.0;
     1268    psF32 midpoint = 0.0;
     1269    psF32 oldMidpoint = 1.0;
     1270    psF32 f = 0.0;
    12711271
    12721272    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
     
    13091309tests to ensure that binNum is within acceptable ranges for both vectors.
    13101310*****************************************************************************/
    1311 float fitQuadraticSearchForYThenReturnX(psVector *xVec,
     1311psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
    13121312                                        psVector *yVec,
    13131313                                        psS32 binNum,
    1314                                         float yVal)
     1314                                        psF32 yVal)
    13151315{
    13161316    PS_VECTOR_CHECK_NULL(xVec, NAN);
     
    13311331    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
    13321332
    1333     float tmpFloat;
     1333    psF32 tmpFloat;
    13341334
    13351335    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
    13361336        // The general case.  We have all three points.
    1337         x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
    1338         x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
    1339         x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
     1337        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
     1338        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
     1339        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
    13401340        y->data.F64[0] = yVec->data.F32[binNum - 1];
    13411341        y->data.F64[1] = yVec->data.F32[binNum];
     
    13461346            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    13471347                    PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
    1348                     (double)yVal,y->data.F64[2],y->data.F64[0]);
     1348                    (psF64)yVal,y->data.F64[2],y->data.F64[0]);
    13491349        }
    13501350        yErr->data.F64[0] = 1.0;
     
    14251425XXX: Check for errors in psLib routines that we call.
    14261426*****************************************************************************/
    1427 int p_psVectorRobustStats(const psVector* restrict myVector,
    1428                           const psVector* restrict errors,
    1429                           const psVector* restrict maskVector,
    1430                           psU32 maskVal,
    1431                           psStats* stats)
     1427psS32 p_psVectorRobustStats(const psVector* restrict myVector,
     1428                            const psVector* restrict errors,
     1429                            const psVector* restrict maskVector,
     1430                            psU32 maskVal,
     1431                            psStats* stats)
    14321432{
    14331433    psHistogram* robustHistogram = NULL;
    14341434    psVector* robustHistogramVector = NULL;
    1435     float binSize = 0.0;        // Size of the histogram bins
     1435    psF32 binSize = 0.0;        // Size of the histogram bins
    14361436    psS32 LQBinNum = -1;          // Bin num for lower quartile
    14371437    psS32 UQBinNum = -1;          // Bin num for upper quartile
     
    14391439    psS32 i = 0;                  // Loop index variable
    14401440    psS32 modeBinNum = 0;
    1441     float modeBinCount = 0.0;
    1442     float dL = 0.0;
     1441    psF32 modeBinCount = 0.0;
     1442    psF32 dL = 0.0;
    14431443    psS32 numBins = 0;
    1444     float myMean = 0.0;
    1445     float myStdev = 0.0;
    1446     float countFloat = 0.0;
    1447     float diff = 0.0;
    1448     float sumSquares = 0.0;
    1449     float sumDiffs = 0.0;
     1444    psF32 myMean = 0.0;
     1445    psF32 myStdev = 0.0;
     1446    psF32 countFloat = 0.0;
     1447    psF32 diff = 0.0;
     1448    psF32 sumSquares = 0.0;
     1449    psF32 sumDiffs = 0.0;
    14501450    psVector* cumulativeRobustSums = NULL;
    1451     float sumRobust = 0.0;
    1452     float sumN50 = 0.0;
    1453     float sumNfit = 0.0;
     1451    psF32 sumRobust = 0.0;
     1452    psF32 sumN50 = 0.0;
     1453    psF32 sumNfit = 0.0;
    14541454    psScalar tmpScalar;
    14551455    tmpScalar.type.type = PS_TYPE_F32;
     
    14601460    // that by 10.0;
    14611461    //XXX: add errors
    1462     int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);
     1462    psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);
    14631463    if (rc != 0) {
    14641464        psError(PS_ERR_UNEXPECTED_NULL,
     
    16141614    for (i=0;i<robustHistogramVector->n;i++) {
    16151615        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1616         ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
     1616        ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    16171617        y->data.F32[i] = robustHistogramVector->data.F32[i];
    16181618    }
     
    17701770    The histogram structure
    17711771 *****************************************************************************/
    1772 psHistogram* psHistogramAlloc(float lower, float upper, psS32 n)
     1772psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n)
    17731773{
    17741774    PS_INT_CHECK_POSITIVE(n, NULL);
     
    17771777    psS32 i = 0;                  // Loop index variable
    17781778    psHistogram* newHist = NULL;        // The new histogram structure
    1779     float binSize = 0.0;        // The histogram bin size
     1779    psF32 binSize = 0.0;        // The histogram bin size
    17801780
    17811781    // Allocate memory for the new histogram structure.  If there are N
     
    17871787
    17881788    // Calculate the bounds for each bin.
    1789     binSize = (upper - lower) / (float)n;
     1789    binSize = (upper - lower) / (psF32)n;
    17901790    // XXX: Is the following necessary? It prevents the max data point
    17911791    // from being in a non-existant bin.
    17921792    binSize += FLT_EPSILON;
    17931793    for (i = 0; i < n + 1; i++) {
    1794         newHist->bounds->data.F32[i] = lower + (binSize * (float)i);
     1794        newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i);
    17951795    }
    17961796
     
    18571857}
    18581858
     1859/*****************************************************************************
     1860UpdateHistogramBins(binNum, out, data, error): This routine is to be used when
     1861updating the histogram in the presence of errors in the input data.  We treat
     1862the data point as a boxcar PDF and update a range of points surrounding the
     1863histogram bin which contains the point.  The width of that boxcar is defined
     1864as 2.35 * error.  Inputs:
     1865    binNum: the bin number of the data point in the histogram
     1866    out: the histogram structure
     1867    data: the data point value
     1868    error: the error in that data point
     1869 
     1870XXX: Must test this.
     1871 *****************************************************************************/
    18591872psS32 UpdateHistogramBins(psS32 binNum,
    18601873                          psHistogram* out,
     
    18631876{
    18641877    PS_PTR_CHECK_NULL(out, -1);
    1865     PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2);
    1866     /*
    1867         psF32 width = 2.35 * error;
    1868         psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum]
    1869         psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0;
    1870      
    1871         if (width <= centerBinWidth) {
    1872             out->nums->data.F32[binNum]+= 1.0;
    1873         } else {
    1874             out->nums->data.F32[binNum]+= centerBinWidth / width;
    1875             // XXX: walk to the left, adding fractional values.
    1876             // XXX: walk to the right, adding fractional values.
    1877      
    1878      
    1879         }
    1880     */
     1878    PS_PTR_CHECK_NULL(out->bounds, -1);
     1879    PS_PTR_CHECK_NULL(out->nums, -1);
     1880    PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2);
     1881    PS_FLOAT_COMPARE(0.0, error, -3);
     1882    PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);
     1883
     1884    psF32 boxcarWidth = 2.35 * error;
     1885    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
     1886                          out->bounds->data.F32[binNum+1]) / 2.0;
     1887    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0);
     1888    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0);
     1889    psS32 bin;
     1890    psS32 boxcarLeftBinNum;
     1891    psS32 boxcarRightBinNum;
     1892
     1893    // Determine the left endpoint of the boxcar for the PDF.
     1894    for (bin=binNum ; bin >= 0 ; bin--) {
     1895        if (out->nums->data.F32[bin] <= boxcarLeft) {
     1896            boxcarLeftBinNum = bin;
     1897            break;
     1898        }
     1899    }
     1900
     1901    // Determine the right endpoint of the boxcar for the PDF.
     1902    for (bin=binNum ; bin < out->nums->n ; bin++) {
     1903        if (out->nums->data.F32[bin] >= boxcarRight) {
     1904            boxcarRightBinNum = bin;
     1905            break;
     1906        }
     1907    }
     1908
     1909    //
     1910    // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
     1911    // bin and return.
     1912    //
     1913    if (boxcarLeftBinNum == boxcarRightBinNum) {
     1914        out->nums->data.F32[binNum]+= 1.0;
     1915        return(0);
     1916    }
     1917
     1918    //
     1919    // If we get here, multiple bins must be updated.  We handle the left
     1920    // endpoint, and right endpoint differently.
     1921    //
     1922    out->nums->data.F32[boxcarLeftBinNum]+=
     1923        (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;
     1924
     1925    //
     1926    // Loop through the center bins, if any.
     1927    //
     1928    for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) {
     1929        out->nums->data.F32[bin]+=
     1930            (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
     1931    }
     1932
     1933    //
     1934    // Handle the right endpoint differently.
     1935    //
     1936    out->nums->data.F32[boxcarRightBinNum]+=
     1937        (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;
     1938
     1939    //
     1940    // Return 0 on success.
     1941    //
    18811942    return(0);
    18821943}
     
    19061967{
    19071968    PS_PTR_CHECK_NULL(out, NULL);
     1969    PS_VECTOR_CHECK_NULL(out->bounds, NULL);
    19081970    PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL);
     1971    PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL);
     1972    PS_VECTOR_CHECK_NULL(out->nums, NULL);
    19091973    PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL);
    19101974    PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL);
     
    19201984
    19211985    psS32 i = 0;                  // Loop index variable
    1922     float binSize = 0.0;          // Histogram bin size
     1986    psF32 binSize = 0.0;          // Histogram bin size
    19231987    psS32 binNum = 0;             // A temporary bin number
    19241988    psS32 numBins = 0;            // The total number of bins
    1925     psS32 tmp = 0;
    19261989    psScalar tmpScalar;
    19271990    tmpScalar.type.type = PS_TYPE_F32;
    1928     psVector* inF32;
    1929     psS32 mustFreeTmp = 1;
     1991    psVector* inF32 = NULL;
     1992    psVector* errorsF32 = NULL;
     1993    psS32 mustFreeVectorIn = 1;
     1994    psS32 mustFreeVectorErrors = 1;
     1995
     1996    // Convert input and errors vectors to F32 if necessary.
    19301997    inF32 = p_psConvertToF32((psVector *) in);
    19311998    if (inF32 == NULL) {
    19321999        inF32 = (psVector *) in;
    1933         mustFreeTmp = 0;
     2000        mustFreeVectorIn = 0;
     2001    }
     2002    errorsF32 = p_psConvertToF32((psVector *) errors);
     2003    if (errorsF32 == NULL) {
     2004        errorsF32 = (psVector *) errors;
     2005        mustFreeVectorErrors = 0;
    19342006    }
    19352007
     
    19502022                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
    19512023                    binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
    1952                     if (errors != NULL) {
     2024                    if (errorsF32 != NULL) {
    19532025                        // XXX: Check return codes.
    19542026                        UpdateHistogramBins(binNum, out,
    19552027                                            inF32->data.F32[i],
    1956                                             errors->data.F32[i]);
     2028                                            errorsF32->data.F32[i]);
    19572029                    } else {
    19582030                        // XXX: This if-statement really shouldn't be necessary.
     
    19692041                    // correct bin number requires a bit more work.
    19702042                    tmpScalar.data.F32 = inF32->data.F32[i];
    1971                     tmp = p_psVectorBinDisect(out->bounds, &tmpScalar);
    1972                     if (tmp < 0) {
     2043                    binNum = p_psVectorBinDisect(out->bounds, &tmpScalar);
     2044                    if (binNum < 0) {
    19732045                        psLogMsg(__func__, PS_LOG_WARN,
    19742046                                 "WARNING: psVectorHistogram(): element outside histogram bounds.\n");
    19752047                    } else {
    1976                         if (errors != NULL) {
     2048                        if (errorsF32 != NULL) {
    19772049                            // XXX: Check return codes.
    1978                             UpdateHistogramBins(tmp, out,
     2050                            UpdateHistogramBins(binNum, out,
    19792051                                                inF32->data.F32[i],
    19802052                                                errors->data.F32[i]);
    19812053                        } else {
    1982                             (out->nums->data.F32[tmp])+= 1.0;
     2054                            (out->nums->data.F32[binNum])+= 1.0;
    19832055                        }
    19842056                    }
     
    19882060    }
    19892061
    1990     if (mustFreeTmp == 1) {
     2062    if (mustFreeVectorIn == 1) {
    19912063        psFree(inF32);
     2064    }
     2065    if (mustFreeVectorErrors == 1) {
     2066        psFree(errorsF32);
    19922067    }
    19932068    return (out);
     
    20152090        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20162091        for (i = 0; i < in->n; i++) {
    2017             tmp->data.F32[i] = (float)in->data.S8[i];
     2092            tmp->data.F32[i] = (psF32)in->data.S8[i];
     2093        }
     2094    } else if (in->type.type == PS_TYPE_S16) {
     2095        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2096        for (i = 0; i < in->n; i++) {
     2097            tmp->data.F32[i] = (psF32) in->data.S16[i];
     2098        }
     2099    } else if (in->type.type == PS_TYPE_S32) {
     2100        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2101        for (i = 0; i < in->n; i++) {
     2102            tmp->data.F32[i] = (psF32)in->data.S32[i];
     2103        }
     2104    } else if (in->type.type == PS_TYPE_S64) {
     2105        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2106        for (i = 0; i < in->n; i++) {
     2107            tmp->data.F32[i] = (psF32)in->data.S64[i];
     2108        }
     2109    } else if (in->type.type == PS_TYPE_U8) {
     2110        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2111        for (i = 0; i < in->n; i++) {
     2112            tmp->data.F32[i] = (psF32)in->data.U8[i];
    20182113        }
    20192114    } else if (in->type.type == PS_TYPE_U16) {
    20202115        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20212116        for (i = 0; i < in->n; i++) {
    2022             tmp->data.F32[i] = (float)in->data.U16[i];
    2023         }
    2024     } else if (in->type.type == PS_TYPE_U8) {
     2117            tmp->data.F32[i] = (psF32)in->data.U16[i];
     2118        }
     2119    } else if (in->type.type == PS_TYPE_U32) {
    20252120        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20262121        for (i = 0; i < in->n; i++) {
    2027             tmp->data.F32[i] = (float)in->data.U8[i];
     2122            tmp->data.F32[i] = (psF32)in->data.U32[i];
     2123        }
     2124    } else if (in->type.type == PS_TYPE_U64) {
     2125        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
     2126        for (i = 0; i < in->n; i++) {
     2127            tmp->data.F32[i] = (psF32)in->data.U64[i];
    20282128        }
    20292129    } else if (in->type.type == PS_TYPE_F64) {
    20302130        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
    20312131        for (i = 0; i < in->n; i++) {
    2032             tmp->data.F32[i] = (float)in->data.F64[i];
     2132            tmp->data.F32[i] = (psF32)in->data.F64[i];
    20332133        }
    20342134    } else if (in->type.type == PS_TYPE_F32) {
    20352135        // do nothing
    20362136    } else {
    2037         char* strType;
    2038         PS_TYPE_NAME(strType,in->type.type);
     2137        psS8* strType;
     2138        PS_TYPE_NAME(strType, in->type.type);
    20392139        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    20402140                PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED,
     
    20742174    }
    20752175
    2076     psVector* inF32;
    2077     psVector* errorsF32;
     2176    psVector* inF32 = NULL;
     2177    psVector* errorsF32 = NULL;
    20782178    psS32 mustFreeVectorIn = 1;
    20792179    psS32 mustFreeVectorErrors = 1;
     
    21782278    return (stats);
    21792279}
     2280
  • trunk/psLib/src/math/psStats.h

    r2778 r2788  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-12-21 20:42:07 $
     12 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-12-22 05:09:32 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5656typedef struct
    5757{
    58     double sampleMean;          ///< formal mean of sample
    59     double sampleMedian;        ///< formal median of sample
    60     double sampleStdev;         ///< standard deviation of sample
    61     double sampleUQ;            ///< upper quartile of sample
    62     double sampleLQ;            ///< lower quartile of sample
    63     double robustMean;          ///< robust mean of array
    64     double robustMedian;        ///< robust median of array
    65     double robustMode;          ///< Robust mode of array
    66     double robustStdev;         ///< robust standard deviation of array
    67     double robustUQ;            ///< robust upper quartile
    68     double robustLQ;            ///< robust lower quartile
     58    psF64 sampleMean;          ///< formal mean of sample
     59    psF64 sampleMedian;        ///< formal median of sample
     60    psF64 sampleStdev;         ///< standard deviation of sample
     61    psF64 sampleUQ;            ///< upper quartile of sample
     62    psF64 sampleLQ;            ///< lower quartile of sample
     63    psF64 robustMean;          ///< robust mean of array
     64    psF64 robustMedian;        ///< robust median of array
     65    psF64 robustMode;          ///< Robust mode of array
     66    psF64 robustStdev;         ///< robust standard deviation of array
     67    psF64 robustUQ;            ///< robust upper quartile
     68    psF64 robustLQ;            ///< robust lower quartile
    6969    psS32 robustN50;              ///<
    7070    psS32 robustNfit;             ///<
    71     double clippedMean;         ///< Nsigma clipped mean
    72     double clippedStdev;        ///< standard deviation after clipping
     71    psF64 clippedMean;         ///< Nsigma clipped mean
     72    psF64 clippedStdev;        ///< standard deviation after clipping
    7373    psS32 clippedNvalues;         ///< ???
    74     double clipSigma;           ///< Nsigma used for clipping; user input
     74    psF64 clipSigma;           ///< Nsigma used for clipping; user input
    7575    psS32 clipIter;               ///< Number of clipping iterations; user input
    76     double min;                 ///< minimum data value in array
    77     double max;                 ///< maximum data value in array
    78     double binsize;             ///<
     76    psF64 min;                 ///< minimum data value in array
     77    psF64 max;                 ///< maximum data value in array
     78    psF64 binsize;             ///<
    7979    psStatsOptions options;     ///< bitmask of calculated values
    8080}
     
    132132 */
    133133psHistogram* psHistogramAlloc(
    134     float lower,                       ///< Lower limit for the bins
    135     float upper,                       ///< Upper limit for the bins
     134    psF32 lower,                       ///< Lower limit for the bins
     135    psF32 upper,                       ///< Upper limit for the bins
    136136    psS32 n                              ///< Number of bins
    137137);
     
    172172    ///< the statistic struct to operate on
    173173
    174     double *value
     174    psF64 *value
    175175    ///< if return is true, this is set to the specified statistic value by stats->options
    176176);
Note: See TracChangeset for help on using the changeset viewer.