Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 1292) +++ /trunk/psLib/src/astro/psCoord.c (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,5 +22,7 @@ #include "psCoord.h" #include "psMemory.h" - +#include "psAbort.h" +#include +#include psPlane *psPlaneTransformApply(psPlane *out, @@ -127,4 +129,117 @@ } +// XXX: Is this the correct way to calculate this? +float cot(float x) +{ + return(1.0 / atan(x)); +} + +float arg(float x, float y) +{ + if (x > 0) { + return(atan(y/x)); + } else if ((x==0) && (y==0)) { + return(0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(-0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(M_PI + atan(y/x)); + } else if ((x==0) && (y==0)) { + return(-M_PI + atan(y/x)); + } + + psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y); + return(0.0); +} + + +psPlane *psProject(const psSphere *coord, + const psProjection *projection) +{ + float R = 0.0; + float alpha = 0.0; + psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane)); + + if (projection->type == PS_PROJ_TAN) { + R = cot(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_SIN) { + R = cos(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->x = coord->d; + tmp->y = coord->r; + + } else if (projection->type == PS_PROJ_MER) { + tmp->x = coord->d; + tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI; + + } else if (projection->type == PS_PROJ_AIT) { + alpha = 1.0 / ((180.0 / M_PI) * + sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5))); + + tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d); + tmp->y = alpha * sin(coord->d); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} + +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection) +{ + float R = 0.0; + float chu = 0.0; + float chu1 = 0.0; + float chu2 = 0.0; + psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere)); + + if (projection->type == PS_PROJ_TAN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = atan(180.0 / (R * M_PI)); + + } else if (projection->type == PS_PROJ_SIN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = acos((R * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->d = coord->x; + tmp->r = coord->y; + + } else if (projection->type == PS_PROJ_MER) { + tmp->d = coord->x; + tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0; + + } else if (projection->type == PS_PROJ_AIT) { + chu1 = (coord->x * M_PI) / 720.0; + chu1*= chu1; + chu2 = (coord->y * M_PI) / 360.0; + chu2*= chu2; + chu = sqrt(1.0 - chu1 - chu2); + tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, + (coord->x * chu * M_PI)/360.0); + tmp->r = asin((coord->y * chu * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} psSphere *psSphereGetOffset(const psSphere *restrict position1, @@ -133,8 +248,11 @@ psSphereOffsetUnit unit) { - return(NULL); -} - - + psAbort(__func__, "This function has not be dedfined.\n"); + return(NULL); +} + + +// XXX: I copied the algorithm from the ADD exactly. Arguments mode and unit +// are ignored. psSphere *psSphereSetOffset(const psSphere *restrict position, const psSphere *restrict offset, @@ -142,4 +260,24 @@ psSphereOffsetUnit unit) { - return(NULL); -} + psPlane lin; + psSphere *tmp; + psProjection proj; + + if (mode == PS_LINEAR) { + proj.R = position->r; + proj.D = position->d; + proj.Xs = 0.0; + proj.Ys = 0.0; + proj.type = PS_PROJ_TAN; + + lin.x = offset->r; + lin.y = offset->d; + + tmp = psDeproject(&lin, &proj); + return(tmp); + } else if (mode == PS_SPHERICAL) { + psAbort(__func__, "Sperical offset modes are not defined.\n"); + } + psAbort(__func__, "Unrecognized offset mode\n"); + return(NULL); +} Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 1292) +++ /trunk/psLib/src/astro/psCoord.h (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,6 +27,6 @@ typedef struct { - double x; ///< x position - double y; ///< y position + double x; ///< x position + double y; ///< y position double xErr; ///< Error in x position double yErr; ///< Error in y position @@ -36,6 +36,6 @@ typedef struct { - double r; ///< RA - double d; ///< Dec + double r; ///< RA + double d; ///< Dec double rErr; ///< Error in RA double dErr; ///< Error in Dec @@ -63,22 +63,24 @@ double sinNPlat; ///< sin of North Pole latitude double cosNPlat; ///< cos of North Pole latitude - double sinZP; ///< sin of Forst PT os Ares lon - double cosZP; ///< cos of Forst PT os Ares lon + double sinZP; ///< sin of Forst PT os Ares lon + double cosZP; ///< cos of Forst PT os Ares lon } psSphereTransform; typedef enum { - PS_PROJ_TAN, ///< Tangent projection - PS_PROJ_SIN, ///< Sine projection - PS_PROJ_AIT, ///< Aitoff projection - PS_PROJ_PAR, ///< Par projection - PS_PROJ_GLS, ///< GLS projection - PS_PROJ_NTYPE ///< Number of types; must be last. + PS_PROJ_TAN, ///< Tangent projection + PS_PROJ_SIN, ///< Sine projection + PS_PROJ_AIT, ///< Aitoff projection + PS_PROJ_PAR, ///< Par projection + PS_PROJ_GLS, ///< GLS projection + PS_PROJ_CAR, ///< CAR projection + PS_PROJ_MER, ///< MER projection + PS_PROJ_NTYPE ///< Number of types; must be last. } psProjectionType; typedef struct { - double R; ///< Coordinates of projection center - double D; ///< Coordinates of projection center + double R; ///< Coordinates of projection center + double D; ///< Coordinates of projection center double Xs; ///< plate-scale in X direction double Ys; ///< plate-scale in Y direction @@ -128,6 +130,6 @@ const psProjection *projection); -psSphere *psProject(const psPlane *coord, - const psProjection *projection); +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection); psSphere *psSphereGetOffset(const psSphere *restrict position1, Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 1292) +++ /trunk/psLib/src/astronomy/psCoord.c (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,5 +22,7 @@ #include "psCoord.h" #include "psMemory.h" - +#include "psAbort.h" +#include +#include psPlane *psPlaneTransformApply(psPlane *out, @@ -127,4 +129,117 @@ } +// XXX: Is this the correct way to calculate this? +float cot(float x) +{ + return(1.0 / atan(x)); +} + +float arg(float x, float y) +{ + if (x > 0) { + return(atan(y/x)); + } else if ((x==0) && (y==0)) { + return(0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(-0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(M_PI + atan(y/x)); + } else if ((x==0) && (y==0)) { + return(-M_PI + atan(y/x)); + } + + psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y); + return(0.0); +} + + +psPlane *psProject(const psSphere *coord, + const psProjection *projection) +{ + float R = 0.0; + float alpha = 0.0; + psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane)); + + if (projection->type == PS_PROJ_TAN) { + R = cot(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_SIN) { + R = cos(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->x = coord->d; + tmp->y = coord->r; + + } else if (projection->type == PS_PROJ_MER) { + tmp->x = coord->d; + tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI; + + } else if (projection->type == PS_PROJ_AIT) { + alpha = 1.0 / ((180.0 / M_PI) * + sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5))); + + tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d); + tmp->y = alpha * sin(coord->d); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} + +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection) +{ + float R = 0.0; + float chu = 0.0; + float chu1 = 0.0; + float chu2 = 0.0; + psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere)); + + if (projection->type == PS_PROJ_TAN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = atan(180.0 / (R * M_PI)); + + } else if (projection->type == PS_PROJ_SIN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = acos((R * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->d = coord->x; + tmp->r = coord->y; + + } else if (projection->type == PS_PROJ_MER) { + tmp->d = coord->x; + tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0; + + } else if (projection->type == PS_PROJ_AIT) { + chu1 = (coord->x * M_PI) / 720.0; + chu1*= chu1; + chu2 = (coord->y * M_PI) / 360.0; + chu2*= chu2; + chu = sqrt(1.0 - chu1 - chu2); + tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, + (coord->x * chu * M_PI)/360.0); + tmp->r = asin((coord->y * chu * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} psSphere *psSphereGetOffset(const psSphere *restrict position1, @@ -133,8 +248,11 @@ psSphereOffsetUnit unit) { - return(NULL); -} - - + psAbort(__func__, "This function has not be dedfined.\n"); + return(NULL); +} + + +// XXX: I copied the algorithm from the ADD exactly. Arguments mode and unit +// are ignored. psSphere *psSphereSetOffset(const psSphere *restrict position, const psSphere *restrict offset, @@ -142,4 +260,24 @@ psSphereOffsetUnit unit) { - return(NULL); -} + psPlane lin; + psSphere *tmp; + psProjection proj; + + if (mode == PS_LINEAR) { + proj.R = position->r; + proj.D = position->d; + proj.Xs = 0.0; + proj.Ys = 0.0; + proj.type = PS_PROJ_TAN; + + lin.x = offset->r; + lin.y = offset->d; + + tmp = psDeproject(&lin, &proj); + return(tmp); + } else if (mode == PS_SPHERICAL) { + psAbort(__func__, "Sperical offset modes are not defined.\n"); + } + psAbort(__func__, "Unrecognized offset mode\n"); + return(NULL); +} Index: /trunk/psLib/src/astronomy/psCoord.h =================================================================== --- /trunk/psLib/src/astronomy/psCoord.h (revision 1292) +++ /trunk/psLib/src/astronomy/psCoord.h (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,6 +27,6 @@ typedef struct { - double x; ///< x position - double y; ///< y position + double x; ///< x position + double y; ///< y position double xErr; ///< Error in x position double yErr; ///< Error in y position @@ -36,6 +36,6 @@ typedef struct { - double r; ///< RA - double d; ///< Dec + double r; ///< RA + double d; ///< Dec double rErr; ///< Error in RA double dErr; ///< Error in Dec @@ -63,22 +63,24 @@ double sinNPlat; ///< sin of North Pole latitude double cosNPlat; ///< cos of North Pole latitude - double sinZP; ///< sin of Forst PT os Ares lon - double cosZP; ///< cos of Forst PT os Ares lon + double sinZP; ///< sin of Forst PT os Ares lon + double cosZP; ///< cos of Forst PT os Ares lon } psSphereTransform; typedef enum { - PS_PROJ_TAN, ///< Tangent projection - PS_PROJ_SIN, ///< Sine projection - PS_PROJ_AIT, ///< Aitoff projection - PS_PROJ_PAR, ///< Par projection - PS_PROJ_GLS, ///< GLS projection - PS_PROJ_NTYPE ///< Number of types; must be last. + PS_PROJ_TAN, ///< Tangent projection + PS_PROJ_SIN, ///< Sine projection + PS_PROJ_AIT, ///< Aitoff projection + PS_PROJ_PAR, ///< Par projection + PS_PROJ_GLS, ///< GLS projection + PS_PROJ_CAR, ///< CAR projection + PS_PROJ_MER, ///< MER projection + PS_PROJ_NTYPE ///< Number of types; must be last. } psProjectionType; typedef struct { - double R; ///< Coordinates of projection center - double D; ///< Coordinates of projection center + double R; ///< Coordinates of projection center + double D; ///< Coordinates of projection center double Xs; ///< plate-scale in X direction double Ys; ///< plate-scale in Y direction @@ -128,6 +130,6 @@ const psProjection *projection); -psSphere *psProject(const psPlane *coord, - const psProjection *projection); +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection); psSphere *psSphereGetOffset(const psSphere *restrict position1, Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 1292) +++ /trunk/psLib/src/dataManip/psStats.c (revision 1293) @@ -131,4 +131,8 @@ mean of the input vector. Inputs + myVector + maskVector + maskVal + stats Returns NULL @@ -404,28 +408,30 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - psStats *stats2 = NULL; // Temporary stats structure // Determine if the number of data points exceed a threshold which will // cause to generate robust stats, as opposed to exact stats. - - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } + // XXX: This code is no longer used. We now calculate the median exactly + // regardless of the vector size. + /* + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleMedian = stats2->robustMedian; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + */ // Determine how many data points fit inside this min/max range @@ -483,4 +489,5 @@ // Calculate the median exactly. + // XXX: Is this the correct action? if (0 == (nValues % 2)) { stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + @@ -498,4 +505,12 @@ This routine smoothes the data in the input robustHistogram with a Gaussian of width sigma. + + XXX: Must consult with IfA on the proper values for the Gaussian + smoothing. Currently, the Gaussian coefficients are such that only + the center coefficient is non-zero. This is because "e" is being + raised to a very large negative number for all points except the + center point. + + XXX: use a static variable for gaussianCoefs[] and compute them once. *****************************************************************************/ psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, @@ -541,4 +556,5 @@ } + // Perform the actual smoothing. for(i=0;inums->n;i++) { smooth->data.F32[i] = 0.0; @@ -576,26 +592,4 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - psStats *stats2 = NULL; - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleUQ = stats2->robustUQ; - stats->sampleLQ = stats2->robustLQ; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } // Determine how many data points fit inside this min/max range @@ -861,4 +855,6 @@ /***************************************************************************** +p_psNormalizeVector(myData): this is a private function which normalizes the +elements of a vector to a range between 0.0 and 1.0. *****************************************************************************/ void p_psNormalizeVector(psVector *myData) @@ -885,4 +881,9 @@ +/***************************************************************************** +p_psGaussian(myData, myParams): an internal function, used by robust stats, +intended to compute the Gaussian with a specified sigma/mean, at the +specified data point. + *****************************************************************************/ float p_psGaussian(const psVector *restrict myData, const psVector *restrict myParams) @@ -894,8 +895,12 @@ tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); + // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); return(tmp); } +/***************************************************************************** +p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above Gaussian function. + *****************************************************************************/ float p_psGaussianDeriv(const psVector *restrict myData, const psVector *restrict myParams, @@ -922,4 +927,9 @@ +/***************************************************************************** +p_psQuadratic(myData, myParams): an internal function, used by robust stats, +intended to compute a quadratic, with the specified parameters, at the +specified data point. + *****************************************************************************/ float p_psQuadratic(const psVector *restrict myParams, const psVector *restrict myCoords) @@ -935,4 +945,8 @@ } +/***************************************************************************** +p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above quadratic function. + *****************************************************************************/ float p_psQuadraticDeriv(const psVector *restrict myParams, const psVector *restrict myCoords, @@ -959,6 +973,6 @@ [rangeLow, rangeHigh]. It determines the x-value of that polynomial such that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. +range and assumes that the polynomial is monotonically increasing or +decreasing within that range. *****************************************************************************/ float p_ps1DPolyMedian(psPolynomial1D *myPoly, @@ -998,49 +1012,58 @@ /****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes -as input a 1-D polynomial of arbitrary order (though we are using 2nd-order -polynomials here) and a range of x-values for which it is defined: -[rangeLow, rangeHigh]. It determines the x-value of that polynomial such -that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram *robustHistogram, +p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given +the specified histogram, with the specified pre-calculated cumulativeSums, +this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", +and then finds the point x such that f(x) == fitFloat. + +XXX: This function is currently not being used. + *****************************************************************************/ +float p_psFitQuadratic(psHistogram *histogram, + psVector *cumulativeSums, int binNum, float fitFloat) { - /* - if ((binNum > 0) && - (binNum < (robustHistogram->nums->n+1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum-1] + - robustHistogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum] + - robustHistogram->bounds->data.F32[binNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum+2]); - - y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1]; - y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum]; - y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1]; - - if (!((y->data.F64[0] <= fitFloat) && - (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat); - } else { - return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum])); - } - */ + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + + if ((binNum > 0) && + (binNum < (histogram->nums->n+1))) { + x->data.F64[0] = (double) 0.5 * + (histogram->bounds->data.F32[binNum-1] + + histogram->bounds->data.F32[binNum]); + x->data.F64[1] = (double) 0.5 * + (histogram->bounds->data.F32[binNum] + + histogram->bounds->data.F32[binNum+1]); + x->data.F64[2] = (double) 0.5 * + (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum+2]); + + y->data.F64[0] = cumulativeSums->data.F32[binNum-1]; + y->data.F64[1] = cumulativeSums->data.F32[binNum]; + y->data.F64[2] = cumulativeSums->data.F32[binNum+1]; + + if (!((y->data.F64[0] <= fitFloat) && + (fitFloat <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], + fitFloat)); + } else { + return(0.5 * (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum])); + } + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myPoly); return(0.0); } @@ -1049,25 +1072,24 @@ p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure calculates a variety of robust stat measures: -PS_STAT_ROBUST_MEAN -PS_STAT_ROBUST_MEDIAN -PS_STAT_ROBUST_MODE -PS_STAT_ROBUST_STDEV -PS_STAT_ROBUST_QUARTILE + PS_STAT_ROBUST_MEAN + PS_STAT_ROBUST_MEDIAN + PS_STAT_ROBUST_MODE + PS_STAT_ROBUST_STDEV + PS_STAT_ROBUST_QUARTILE I have included all that computation in a single function, as opposed to -breaking it across several functions for one primary reason: -they all +breaking it across several functions for one primary reason: they all require the same basic initial processing steps (calculate the histogram, -etc.) -however there is no currently defined means, in the SDRS, to -specify this computation as a separate step. If the above robust stat -measures were calcualted in separate functiosn, then much of the initial -processing would be duplicated. +etc.) however there is no currently defined means, in the SDRS, to specify +this computation as a separate step. If the above robust stat measures were +calcualted in separate functiosn, then much of the initial processing would +be duplicated. + Inputs -myVector -maskVector -maskVal -stats + myVector + maskVector + maskVal + stats Returns -NULL + NULL *****************************************************************************/ void p_psVectorRobustStats(const psVector *restrict myVector, @@ -1224,4 +1246,5 @@ sumN50+= (float) robustHistogram->nums->data.S32[i]; } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? dL = (UQBinNum - LQBinNum) / 4; @@ -1244,5 +1267,5 @@ robustHistogramVector->data.F32[i]; cumulativeMedian+= robustHistogramVector->data.F32[i]; - sumNfit+= (float) robustHistogramVector->data.S32[i]; + sumNfit+= (float) robustHistogram->nums->data.S32[i]; } } @@ -1780,4 +1803,12 @@ } + // XXX: Should we abort if (stats->min == stats->max)? + if (stats->min >= stats->max) { + psAbort(__func__, "psVectorStats() called with range: %f to %f\n", + stats->min, stats->max); + } + + + PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); if (mask != NULL) { Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 1292) +++ /trunk/psLib/src/math/psStats.c (revision 1293) @@ -131,4 +131,8 @@ mean of the input vector. Inputs + myVector + maskVector + maskVal + stats Returns NULL @@ -404,28 +408,30 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - psStats *stats2 = NULL; // Temporary stats structure // Determine if the number of data points exceed a threshold which will // cause to generate robust stats, as opposed to exact stats. - - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } + // XXX: This code is no longer used. We now calculate the median exactly + // regardless of the vector size. + /* + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleMedian = stats2->robustMedian; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + */ // Determine how many data points fit inside this min/max range @@ -483,4 +489,5 @@ // Calculate the median exactly. + // XXX: Is this the correct action? if (0 == (nValues % 2)) { stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + @@ -498,4 +505,12 @@ This routine smoothes the data in the input robustHistogram with a Gaussian of width sigma. + + XXX: Must consult with IfA on the proper values for the Gaussian + smoothing. Currently, the Gaussian coefficients are such that only + the center coefficient is non-zero. This is because "e" is being + raised to a very large negative number for all points except the + center point. + + XXX: use a static variable for gaussianCoefs[] and compute them once. *****************************************************************************/ psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, @@ -541,4 +556,5 @@ } + // Perform the actual smoothing. for(i=0;inums->n;i++) { smooth->data.F32[i] = 0.0; @@ -576,26 +592,4 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - psStats *stats2 = NULL; - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleUQ = stats2->robustUQ; - stats->sampleLQ = stats2->robustLQ; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } // Determine how many data points fit inside this min/max range @@ -861,4 +855,6 @@ /***************************************************************************** +p_psNormalizeVector(myData): this is a private function which normalizes the +elements of a vector to a range between 0.0 and 1.0. *****************************************************************************/ void p_psNormalizeVector(psVector *myData) @@ -885,4 +881,9 @@ +/***************************************************************************** +p_psGaussian(myData, myParams): an internal function, used by robust stats, +intended to compute the Gaussian with a specified sigma/mean, at the +specified data point. + *****************************************************************************/ float p_psGaussian(const psVector *restrict myData, const psVector *restrict myParams) @@ -894,8 +895,12 @@ tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); + // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); return(tmp); } +/***************************************************************************** +p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above Gaussian function. + *****************************************************************************/ float p_psGaussianDeriv(const psVector *restrict myData, const psVector *restrict myParams, @@ -922,4 +927,9 @@ +/***************************************************************************** +p_psQuadratic(myData, myParams): an internal function, used by robust stats, +intended to compute a quadratic, with the specified parameters, at the +specified data point. + *****************************************************************************/ float p_psQuadratic(const psVector *restrict myParams, const psVector *restrict myCoords) @@ -935,4 +945,8 @@ } +/***************************************************************************** +p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above quadratic function. + *****************************************************************************/ float p_psQuadraticDeriv(const psVector *restrict myParams, const psVector *restrict myCoords, @@ -959,6 +973,6 @@ [rangeLow, rangeHigh]. It determines the x-value of that polynomial such that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. +range and assumes that the polynomial is monotonically increasing or +decreasing within that range. *****************************************************************************/ float p_ps1DPolyMedian(psPolynomial1D *myPoly, @@ -998,49 +1012,58 @@ /****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes -as input a 1-D polynomial of arbitrary order (though we are using 2nd-order -polynomials here) and a range of x-values for which it is defined: -[rangeLow, rangeHigh]. It determines the x-value of that polynomial such -that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram *robustHistogram, +p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given +the specified histogram, with the specified pre-calculated cumulativeSums, +this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", +and then finds the point x such that f(x) == fitFloat. + +XXX: This function is currently not being used. + *****************************************************************************/ +float p_psFitQuadratic(psHistogram *histogram, + psVector *cumulativeSums, int binNum, float fitFloat) { - /* - if ((binNum > 0) && - (binNum < (robustHistogram->nums->n+1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum-1] + - robustHistogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum] + - robustHistogram->bounds->data.F32[binNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum+2]); - - y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1]; - y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum]; - y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1]; - - if (!((y->data.F64[0] <= fitFloat) && - (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat); - } else { - return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum])); - } - */ + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + + if ((binNum > 0) && + (binNum < (histogram->nums->n+1))) { + x->data.F64[0] = (double) 0.5 * + (histogram->bounds->data.F32[binNum-1] + + histogram->bounds->data.F32[binNum]); + x->data.F64[1] = (double) 0.5 * + (histogram->bounds->data.F32[binNum] + + histogram->bounds->data.F32[binNum+1]); + x->data.F64[2] = (double) 0.5 * + (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum+2]); + + y->data.F64[0] = cumulativeSums->data.F32[binNum-1]; + y->data.F64[1] = cumulativeSums->data.F32[binNum]; + y->data.F64[2] = cumulativeSums->data.F32[binNum+1]; + + if (!((y->data.F64[0] <= fitFloat) && + (fitFloat <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], + fitFloat)); + } else { + return(0.5 * (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum])); + } + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myPoly); return(0.0); } @@ -1049,25 +1072,24 @@ p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure calculates a variety of robust stat measures: -PS_STAT_ROBUST_MEAN -PS_STAT_ROBUST_MEDIAN -PS_STAT_ROBUST_MODE -PS_STAT_ROBUST_STDEV -PS_STAT_ROBUST_QUARTILE + PS_STAT_ROBUST_MEAN + PS_STAT_ROBUST_MEDIAN + PS_STAT_ROBUST_MODE + PS_STAT_ROBUST_STDEV + PS_STAT_ROBUST_QUARTILE I have included all that computation in a single function, as opposed to -breaking it across several functions for one primary reason: -they all +breaking it across several functions for one primary reason: they all require the same basic initial processing steps (calculate the histogram, -etc.) -however there is no currently defined means, in the SDRS, to -specify this computation as a separate step. If the above robust stat -measures were calcualted in separate functiosn, then much of the initial -processing would be duplicated. +etc.) however there is no currently defined means, in the SDRS, to specify +this computation as a separate step. If the above robust stat measures were +calcualted in separate functiosn, then much of the initial processing would +be duplicated. + Inputs -myVector -maskVector -maskVal -stats + myVector + maskVector + maskVal + stats Returns -NULL + NULL *****************************************************************************/ void p_psVectorRobustStats(const psVector *restrict myVector, @@ -1224,4 +1246,5 @@ sumN50+= (float) robustHistogram->nums->data.S32[i]; } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? dL = (UQBinNum - LQBinNum) / 4; @@ -1244,5 +1267,5 @@ robustHistogramVector->data.F32[i]; cumulativeMedian+= robustHistogramVector->data.F32[i]; - sumNfit+= (float) robustHistogramVector->data.S32[i]; + sumNfit+= (float) robustHistogram->nums->data.S32[i]; } } @@ -1780,4 +1803,12 @@ } + // XXX: Should we abort if (stats->min == stats->max)? + if (stats->min >= stats->max) { + psAbort(__func__, "psVectorStats() called with range: %f to %f\n", + stats->min, stats->max); + } + + + PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); if (mask != NULL) {