Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 1784) +++ trunk/psLib/src/math/psPolynomial.c (revision 1823) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 23:20:29 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-17 02:14:52 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,4 +24,5 @@ #include "psMemory.h" #include "psVector.h" +#include "psScalar.h" #include "psTrace.h" #include "psError.h" @@ -1682,14 +1683,16 @@ /***************************************************************************** -VectorBinDisect(): This is a private function which takes as input a vector -of floating point data as well as a single floating point values. The input -vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). -This routine does a binary disection of the vector and returns "i" such -that (v[i] <= x <= v[i+1). If x lies outside the range of v[], +VectorBinDisectF32(): This is a private function which takes as input a +vector of floating point data as well as a single floating point values. +The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for +all j>=i). This routine does a binary disection of the vector and returns +"i" such that (v[i] <= x <= v[i+1). If x lies outside the range of v[], then this routine prints a warning message and returns -1. + +XXX: Macro this for a few different types. *****************************************************************************/ -int VectorBinDisect(float *bins, - int numBins, - float x) +int VectorBinDisectF32(float *bins, + int numBins, + float x) { int min; @@ -1697,11 +1700,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, - "---- Calling VectorBinDisect(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + "---- Calling VectorBinDisectF32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", + "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1713,11 +1716,59 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, - "---- Exiting VectorBinDisect(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + "---- Exiting VectorBinDisectF32(): bin %d\n", mid); + return(mid); + } else if (x < bins[mid]) { + max = mid-1; + } else { + min = mid; + } + mid = ((max+1)+min)/2; + } + + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + "---- Exiting VectorBinDisectF32(): bin %d\n", min); + return(min); +} + +/***************************************************************************** +VectorBinDisectS32(): integer version of above. + *****************************************************************************/ +int VectorBinDisectS32(int *bins, + int numBins, + int x) +{ + int min; + int max; + int mid; + + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + "---- Calling VectorBinDisectS32(%f)\n", x); + + if ((x < bins[0]) || + (x > bins[numBins-1])) { + psLogMsg(__func__, PS_LOG_WARN, + "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", + x, bins[0], bins[numBins-1]); + return(-1); + } + + min = 0; + max = numBins-2; + mid = ((max+1)-min)/2; + + while (min != max) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", + min, mid, max, x, bins[mid]); + + if (x == bins[mid]) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + "---- Exiting VectorBinDisectS32(): bin %d\n", min); return(min); } else if (x < bins[mid]) { @@ -1729,8 +1780,131 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, - "---- Exiting VectorBinDisect(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + "---- Exiting VectorBinDisectS32(): bin %d\n", min); return(min); } + +/***************************************************************************** +p_psVectorBinDisect(): A wrapper to the above VectorBinDisect(). + *****************************************************************************/ +int p_psVectorBinDisect(psVector *bins, + psScalar *x) +{ + if (x->type.type != bins->type.type) { + // XXX: Generate error message. + return(-1); + } + + if (x->type.type == PS_TYPE_S32) { + return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + } else if (x->type.type == PS_TYPE_F32) { + return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + } else { + // XXX: Generate error message. + } + return(-1); +} + +/***************************************************************************** +p_psInterpolate1D(): This routine will take as input n-element floating +point arrays domain and range, and the x value, assumed to lie with the +domain vector. It produces as output the n-order LaGrange interpolated +value of x. + +XXX: do we error check for non-distinct domain values? + *****************************************************************************/ +float p_psFullInterpolate1DF32(float *domain, + float *range, + int n, + float x) +{ + int i; + int m; + static psVector *p = NULL; + psVectorRecycle(p, n, PS_TYPE_F32); + p_psMemSetPersistent(p, true); + + // From NR, during each iteration of the m loop, we are computing the + // p_{i ... i+m} terms. + for (m=0;mdata.F32[i] = range[i]; + } else { + // From NR: we are computing P_{i ... i+m} + p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) + + ((range[i]-x) * p->data.F32[i+1])) / + (domain[i] - domain[i+m]); + } + } + } + return(p->data.F32[0]); +} + + +// This is a +float p_psInterpolate1DF32(float *domain, + float *range, + int n, + int order, + float x) +{ + int binNum; + int numIntPoints = order+1; + int origin; + + binNum = VectorBinDisectF32(domain, n, x); + if (0 == numIntPoints%2) { + origin = binNum - ((numIntPoints/2) - 1); + } else { + origin = binNum - (numIntPoints/2); + if ((x-domain[binNum]) > (domain[binNum+1]-x)) { + // x is closer to binNum+1. + origin = 1 + (binNum - (numIntPoints/2)); + } + } + if (origin < 0) { + origin = 0; + } + if ((origin + numIntPoints) > n) { + origin = n - numIntPoints; + } + + return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x)); +} + +/***************************************************************************** +p_psInterpolate1D(): This routine will take as input psVectors domain and +range, and the x value, assumed to lie with the domain vector. It produces +as output the LaGrange interpolated value of a polynomial of the specified +order around the point x. + +XXX: This stuff does not work with a mask. + *****************************************************************************/ +float p_psInterpolate1D(psVector *domain, + psVector *range, + int order, + psScalar *x) +{ + if (domain->type.type != range->type.type != x->type.type) { + // XXX psError + } + if (domain->n != range->n) { + // XXX psError + } + if (order > (domain->n - 1)) { + // XXX psError: not enough data points for order-order interpolation. + } + + if (x->type.type == PS_TYPE_F32) { + return(p_psInterpolate1DF32(domain->data.F32, range->data.F32, + domain->n, order, x->data.F32)); + } else { + // XXX psError: type not supported + } + + return(-1.0); +} + float psSpline1DEval(const psSpline1D *spline, @@ -1741,5 +1915,5 @@ n = spline->n; - binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); + binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x); if (binNum == -1) { psLogMsg(__func__, PS_LOG_WARN,