样条之Akima光滑插值函数

核心代码：
  1 //////////////////////////////////////////////////////////////////////
  2 // Akima光滑插值
  3 // t    - 存放指定的插值点的值
  4 // s[]  - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，
  5 //        s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）
  6 // k    - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
  7 //////////////////////////////////////////////////////////////////////
  8 static float GetValueAkima(const void* valuesPtr, int stride, int n, float t, float* s, int k)
  9 { 
 10     int kk,m,l;
 11     float u[5],p,q;
 12 
 13     // 初值
 14     memset(s, 0, 5*sizeof(float));
 15 
 16     // 特例处理
 17     if (n < 1) 
 18     {
 19         return s[4];
 20     }
 21     if (n == 1) 
 22     { 
 23         s[4] = YfGetFloatValue(valuesPtr, stride, 0);  
 24         s[0] = s[4];
 25         return s[4];
 26     }
 27 
 28     float xStep = 1.0f/(n - 1);
 29 
 30     if (n == 2)
 31     { 
 32         float y0 = YfGetFloatValue(valuesPtr, stride, 0); 
 33         float y1 = YfGetFloatValue(valuesPtr, stride, 1); 
 34         s[0] = y0;
 35         s[1] = (y1-y0)/xStep;
 36         s[4] = y0 + (y1 - y0)*t;
 37         return s[4];
 38     }
 39 
 40     // 插值
 41     if (k<0)
 42     { 
 43         if (t <= xStep) 
 44             kk = 0;
 45         else if (t >= (n-1)*xStep) 
 46             kk = n-2;
 47         else
 48         { 
 49             kk = 1; 
 50             m = n;
 51             while (((kk-m)!=1)&&((kk-m)!=-1))
 52             { 
 53                 l=(kk+m)/2;
 54                 if (t < (l-1)*xStep) 
 55                     m=l;
 56                 else 
 57                     kk=l;
 58             }
 59 
 60             kk=kk-1;
 61         }
 62     }
 63     else 
 64         kk=k;
 65 
 66     if (kk>=n-1) 
 67         kk=n-2;
 68 
 69     u[2]=(YfGetFloatValue(valuesPtr, stride, kk+1) - YfGetFloatValue(valuesPtr, stride, kk))/xStep;
 70     if (n==3)
 71     { 
 72         if (kk==0)
 73         { 
 74             u[3]=(YfGetFloatValue(valuesPtr, stride, 2) - YfGetFloatValue(valuesPtr, stride, 1))/xStep;
 75             u[4]=2.0f*u[3]-u[2];
 76             u[1]=2.0f*u[2]-u[3];
 77             u[0]=2.0f*u[1]-u[2];
 78         }
 79         else
 80         { 
 81             u[1]=(YfGetFloatValue(valuesPtr, stride, 1) - YfGetFloatValue(valuesPtr, stride, 0))/xStep;
 82             u[0]=2.0f*u[1]-u[2];
 83             u[3]=2.0f*u[2]-u[1];
 84             u[4]=2.0f*u[3]-u[2];
 85         }
 86     }
 87     else
 88     { 
 89         if (kk<=1)
 90         { 
 91             u[3]=(YfGetFloatValue(valuesPtr, stride, kk+2) - YfGetFloatValue(valuesPtr, stride, kk+1))/xStep;
 92             if (kk==1)
 93             { 
 94                 u[1]=(YfGetFloatValue(valuesPtr, stride, 1) - YfGetFloatValue(valuesPtr, stride, 0))/xStep;
 95                 u[0]=2.0f*u[1]-u[2];
 96                 if (n==4) 
 97                     u[4]=2.0f*u[3]-u[2];
 98                 else 
 99                     u[4]=(YfGetFloatValue(valuesPtr, stride, 4) - YfGetFloatValue(valuesPtr, stride, 3))/xStep;
100             }
101             else
102             { 
103                 u[1]=2.0f*u[2]-u[3];
104                 u[0]=2.0f*u[1]-u[2];
105                 u[4]=(YfGetFloatValue(valuesPtr, stride, 3) - YfGetFloatValue(valuesPtr, stride, 2))/xStep;
106             }
107         }
108         else if (kk>=(n-3))
109         { 
110             u[1]=(YfGetFloatValue(valuesPtr, stride, kk) - YfGetFloatValue(valuesPtr, stride, kk-1))/xStep;
111             if (kk==(n-3))
112             { 
113                 u[3]=(YfGetFloatValue(valuesPtr, stride, n-1) - YfGetFloatValue(valuesPtr, stride, n-2))/xStep;
114                 u[4]=2.0f*u[3]-u[2];
115                 if (n==4) 
116                     u[0]=2.0f*u[1]-u[2];
117                 else 
118                     u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
119             }
120             else
121             { 
122                 u[3]=2.0f*u[2]-u[1];
123                 u[4]=2.0f*u[3]-u[2];
124                 u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
125             }
126         }
127         else
128         { 
129             u[1]=(YfGetFloatValue(valuesPtr, stride, kk  ) - YfGetFloatValue(valuesPtr, stride, kk-1))/xStep;
130             u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
131             u[3]=(YfGetFloatValue(valuesPtr, stride, kk+2) - YfGetFloatValue(valuesPtr, stride, kk+1))/xStep;
132             u[4]=(YfGetFloatValue(valuesPtr, stride, kk+3) - YfGetFloatValue(valuesPtr, stride, kk+2))/xStep;
133         }
134     }
135 
136     s[0] = fabs(u[3]-u[2]);
137     s[1] = fabs(u[0]-u[1]);
138     if ((s[0]+1.0f == 1.0f) && (s[1]+1.0f == 1.0f))
139         p = (u[1]+u[2])/2.0f;
140     else 
141         p = (s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
142 
143     s[0] = fabs(u[3]-u[4]);
144     s[1] = fabs(u[2]-u[1]);
145     if ((s[0]+1.0f==1.0f) && (s[1]+1.0f==1.0f))
146         q = (u[2]+u[3])/2.0f;
147     else 
148         q = (s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
149 
150     s[0] = YfGetFloatValue(valuesPtr, stride, kk);
151     s[1] = p;
152     s[3] = xStep;
153     s[2] = (3.0f*u[2]-2.0f*p-q)/s[3];
154     s[3] = (q+p-2.0f*u[2])/(s[3]*s[3]);
155 
156     //if (k<0)
157     { 
158         p=t-(kk*xStep);
159         s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
160     }
161 
162     return s[4];
163 
164 }
切图：
、
相关软件的下载地址为：https://files.cnblogs.com/WhyEngine/TestSpline.zip
posted on 2014-10-18 13:43 叶飞影阅读(3984) 评论(0) 编辑收藏举报