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             else102             { 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             else121             { 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         else128         { 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 }