IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 11, 2004, 9:32:20 AM (22 years ago)
Author:
gusciora
Message:

...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/dataManip/psFunctions.c

    r2329 r2338  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-11-10 23:22:32 $
     9 *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-11-11 19:32:20 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    816816
    817817/*****************************************************************************
    818 p_psInterpolate1D(): This routine will take as input n-element floating
     818fullInterpolate1DF32(): This routine will take as input n-element floating
    819819point arrays domain and range, and the x value, assumed to lie with the
    820820domain vector.  It produces as output the (n-1)-order LaGrange interpolated
     
    823823XXX: do we error check for non-distinct domain values?
    824824 *****************************************************************************/
    825 static float fullInterpolate1DF32(float *domain,
    826                                   float *range,
    827                                   psS32 n,
    828                                   float x)
    829 {
    830     PS_INT_CHECK_NON_NEGATIVE(n, NAN);
    831     PS_PTR_CHECK_NULL(domain, NAN);
    832     PS_PTR_CHECK_NULL(range, NAN);
    833 
    834     psS32 i;
    835     psS32 m;
    836     static psVector *p = NULL;
    837     p = psVectorRecycle(p, n, PS_TYPE_F32);
    838     p_psMemSetPersistent(p, true);
    839     p_psMemSetPersistent(p->data.F32, true);
    840     /*
    841         psVector *p = psVectorAlloc(n, PS_TYPE_F32);
    842         float tmp;
    843     */
    844 
    845     psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4,
    846             "---- fullInterpolate1DF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
    847 
    848     for (i=0;i<n;i++) {
    849         psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
    850                 "domain/range is (%f %f)\n", domain[i], range[i]);
    851     }
    852 
    853     for (i=0;i<n;i++) {
    854         p->data.F32[i] = range[i];
    855         psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
    856                 "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
    857 
    858     }
    859 
    860     // From NR, during each iteration of the m loop, we are computing the
    861     // p_{i ... i+m} terms.
    862     for (m=1;m<n;m++) {
    863         for (i=0;i<n-m;i++) {
    864             // From NR: we are computing P_{i ... i+m}
    865             p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
    866                               ((domain[i]-x) * p->data.F32[i+1])) /
    867                              (domain[i] - domain[i+m]);
    868             //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
    869             psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
    870                     "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
    871         }
    872     }
    873     psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4,
    874             "---- fullInterpolate1DF32() end ----\n");
    875 
    876     /*
    877         tmp = p->data.F32[0];
    878         psFree(p);
    879         return(tmp);
    880     */
    881     return(p->data.F32[0]);
    882 }
    883 
     825#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
     826static float fullInterpolate1D##TYPE(float *domain, \
     827                                     float *range, \
     828                                     psS32 n, \
     829                                     float x) \
     830{ \
     831    \
     832    psS32 i; \
     833    psS32 m; \
     834    static psVector *p = NULL; \
     835    p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \
     836    p_psMemSetPersistent(p, true); \
     837    p_psMemSetPersistent(p->data.TYPE, true); \
     838    \
     839    psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \
     840            "---- fullInterpolate1D##TYPE() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n); \
     841    \
     842    for (i=0;i<n;i++) { \
     843        psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \
     844                "domain/range is (%f %f)\n", domain[i], range[i]); \
     845    } \
     846    \
     847    for (i=0;i<n;i++) { \
     848        p->data.TYPE[i] = range[i]; \
     849        psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \
     850                "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \
     851        \
     852    } \
     853    \
     854    /* From NR, during each iteration of the m loop, we are computing the \
     855       p_{i ... i+m} terms. \
     856    */ \
     857    for (m=1;m<n;m++) { \
     858        for (i=0;i<n-m;i++) { \
     859            /* From NR: we are computing P_{i ... i+m} \
     860             */ \
     861            p->data.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \
     862                               ((domain[i]-x) * p->data.TYPE[i+1])) / \
     863                              (domain[i] - domain[i+m]); \
     864            /*printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.TYPE[i], domain[i], x, p->data.TYPE[i+1], domain[i], domain[i+m]); \
     865             */ \
     866            psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \
     867                    "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \
     868        } \
     869    } \
     870    psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \
     871            "---- fullInterpolate1D##TYPE() end ----\n"); \
     872    \
     873    return(p->data.TYPE[0]); \
     874} \
     875
     876FUNC_MACRO_FULL_INTERPOLATE_1D(U8)
     877FUNC_MACRO_FULL_INTERPOLATE_1D(U16)
     878FUNC_MACRO_FULL_INTERPOLATE_1D(U32)
     879FUNC_MACRO_FULL_INTERPOLATE_1D(U64)
     880FUNC_MACRO_FULL_INTERPOLATE_1D(S8)
     881FUNC_MACRO_FULL_INTERPOLATE_1D(S16)
     882FUNC_MACRO_FULL_INTERPOLATE_1D(S32)
     883FUNC_MACRO_FULL_INTERPOLATE_1D(S64)
     884FUNC_MACRO_FULL_INTERPOLATE_1D(F32)
     885FUNC_MACRO_FULL_INTERPOLATE_1D(F64)
    884886
    885887/*****************************************************************************
Note: See TracChangeset for help on using the changeset viewer.