Changeset 2338 for trunk/psLib/src/math/psSpline.c
- Timestamp:
- Nov 11, 2004, 9:32:20 AM (22 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psSpline.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psSpline.c
r2329 r2338 7 7 * polynomials. It also contains a Gaussian functions. 8 8 * 9 * @version $Revision: 1.6 2$ $Name: not supported by cvs2svn $10 * @date $Date: 2004-11-1 0 23:22:32$9 * @version $Revision: 1.63 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-11-11 19:32:20 $ 11 11 * 12 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 816 816 817 817 /***************************************************************************** 818 p_psInterpolate1D(): This routine will take as input n-element floating818 fullInterpolate1DF32(): This routine will take as input n-element floating 819 819 point arrays domain and range, and the x value, assumed to lie with the 820 820 domain vector. It produces as output the (n-1)-order LaGrange interpolated … … 823 823 XXX: do we error check for non-distinct domain values? 824 824 *****************************************************************************/ 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) \ 826 static 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 876 FUNC_MACRO_FULL_INTERPOLATE_1D(U8) 877 FUNC_MACRO_FULL_INTERPOLATE_1D(U16) 878 FUNC_MACRO_FULL_INTERPOLATE_1D(U32) 879 FUNC_MACRO_FULL_INTERPOLATE_1D(U64) 880 FUNC_MACRO_FULL_INTERPOLATE_1D(S8) 881 FUNC_MACRO_FULL_INTERPOLATE_1D(S16) 882 FUNC_MACRO_FULL_INTERPOLATE_1D(S32) 883 FUNC_MACRO_FULL_INTERPOLATE_1D(S64) 884 FUNC_MACRO_FULL_INTERPOLATE_1D(F32) 885 FUNC_MACRO_FULL_INTERPOLATE_1D(F64) 884 886 885 887 /*****************************************************************************
Note:
See TracChangeset
for help on using the changeset viewer.
