C语言实现matlab的interp2，函数

　　项目要用到matlab中的Vq = interp2(X,Y,V,Xq,Yq)函数，即把一个已知经纬度和对应值的矩阵，插值变换到一个给定经纬度网格中，也就是对给定网格填值，需要用到插值，这里使用双线性内插法。

*(这只是一个初步完成代码，仅供参考)

以下是对应C代码和测试程序：

  1 //************************************
  2 // 函数名称: inter_linear()
  3 // 函数说明：计算两点之间某一给定点的值，（x0,y0）->(x1,y1),已知x(x0<x<x1),求y
  4 // 返 回 值: double
  5 // 参    数: double x0,y0,x1,y1/*in*/    两个点的坐标（x0,y0）（x1,y1）
  6 
  7 // 作    者：WSS
  8 // 作成日期：2019/09/05
  9 
 10 //************************************
 11 double inter_linear(double x0, double y0, double x1, double y1, double x)
 12 {
 13     double a0, a1, y;
 14     a0 = (x - x1) / (x0 - x1);
 15     a1 = (x - x0) / (x1 - x0);
 16     y = a0*y0 + a1*y1;
 17     return y;
 18 }

  1 //************************************
  2 // 函数名称: interp2d()
  3 // 函数说明：二维插值，同matlab的interp2()功能
  4 // 返 回 值: double
  5 // 参    数: x,y分别为长度为m和n的向量（一维数组），z为矩阵（对应的二维数组(m,n)）
  6 //           a,b分别为长度为asize和bsize的向量（一维数组），out_result为矩阵（对应的二维数组(asize,bsize)）
  7 
  8 // 作    者：WSS
  9 // 作成日期：2019/09/05
 10 
 11 //************************************
 12 int interp2d(double *x, double *y, double *z, int m, int n, double *a, double *b, int asize,int bsize,double *out_result)
 13 {
 14     //a,b是待插值的矩阵行和列，遍历a，取出行，遍历b取出列，将插值计算后的值存储在out_result中
 15     double pointa = 0;
 16     double pointb = 0;
 17     double w1, w2,w;
 18     const int nu_val = -9999;
 19     int tempi = 0;
 20     int tempj =0;
 21     for (int i = 0; i < asize; ++i)
 22     {
 23         for (int j = 0; j < bsize; ++j)
 24         {
 25             // 找到网格点位置
 26             pointa = a[i];
 27             pointb = b[j];
 28             //确定pointa pointb所在的位置
 29             for (int p = 0; p < m-1; ++p)
 30             {
 31                 if (pointa < x[0]||pointa>x[m-1])//范围外的值不进行插值处理
 32                 {
 33                     tempi = nu_val;
 34                     break;
 35                 }
 36                 else if (pointa == x[m - 1])
 37                 {
 38                     tempi = m - 1;
 39                     break;
 40                 }
 41                 else if(pointa >= x[p] && pointa < x[p+1])//x升序
 42                 {
 43                     tempi = p;
 44                     break;
 45                 }
 46             }
 47             for (int q = 0; q < n - 1; ++q)
 48             {
 49                 if (pointb< y[0] || pointb>y[n- 1])//范围外的值不进行插值处理
 50                 {
 51                     tempj = nu_val;
 52                     break;
 53                 }
 54                 else if (pointb == y[n - 1])
 55                 {
 56                     tempj = n - 1;
 57                     break;
 58                 }
 59                 else if (pointb >= y[q] && pointb < y[q + 1])//y升序
 60                 {
 61                     tempj = q;
 62                     break;
 63                 }
 64             }
 65             if (tempj == nu_val || tempi == nu_val)
 66             {
 67                 out_result[i*bsize + j] = nu_val;
 68             }
 69             else
 70             {
 71                 //如果四个点，有一个点为-9999，就不进行插值了
 72                 if(z[tempi*n + tempj]<-900||z[tempi*n + tempj + 1]<-900||z[(tempi + 1)*n + tempj]<-900||z[(tempi + 1)*n + tempj + 1]<-900)
 73                 {
 74                     out_result[i*bsize + j] = nu_val;
 75                 }
 76                 else
 77                 {
 78                     //cout << tempi << " " << tempj << endl;
 79                     /**************x方向进行插值*************/
 80                     if (x[tempi] == pointa)
 81                     {
 82                         //取网格节点值
 83                         w1 = z[tempi*n + tempj];
 84                         w2 = z[tempi*n + tempj + 1];
 85                         //y方向进行插值
 86                         if (y[tempj] == pointb)
 87                         {
 88                             w = w1;
 89                         }
 90                         else
 91                         {
 92                             w = inter_linear(y[tempj], w1, y[tempj + 1], w2, pointb);
 93                         }
 94                     }
 95                     else
 96                     {
 97                         //x方向进行插值
 98                         w1 = inter_linear(x[tempi], z[tempi*n + tempj], x[tempi + 1], z[(tempi + 1)*n + tempj], pointa);
 99                         w2 = inter_linear(x[tempi], z[tempi*n + tempj + 1], x[tempi + 1], z[(tempi + 1)*n + tempj + 1], pointa);
100                         /*******y方向进行插值********/
101 
102                         if (y[tempj] == pointb)
103 
104                         {
105 
106                             //取网格节点值
107 
108                             w = w1;
109 
110                         }
111                         else
112                         {
113                             //进行插值(y)
114 
115                             w = inter_linear(y[tempj], w1, y[tempj + 1], w2, pointb);
116                         }
117                     }
118                     out_result[i*bsize + j] = w;
119                 }
120 
121             }
122 
123 
124         }
125     }
126     return 0;
127 }

测试函数：

 1 int main()
 2 {
 3     //插值函数测试
 4     int ret = 0;
 5     double a1[5] = { 10, 20, 30, 40, 50 };
 6     double a2[6] = {1,2,4,6,8,9};
 7     double z[30] ;
 8     for (int i = 0; i < 5; ++i)
 9     {
10         for (int j = 0; j < 6; ++j)
11         {
12             z[i*6+j] = a1[i]+a2[j];
13         }
14     }
15 
16     double b1[10] = { 10, 20, 30, 40, 50, 10, 20, 30, 40, 50 };
17     double b2[18] = { 0, 2, 4, 6, 8, 10, 0, 2, 4, 6, 8, 10, 0, 2, 4, 6, 8, 10 };
18     double zout[180];
19     ret = interp2d(a1, a2, z, 5, 6, b1, b2, 10, 18, zout);
20     for (int i = 0; i < 5; ++i)
21     {
22         for (int j = 0; j < 6; j++)
23         {
24             cout << z[i * 6 + j] << " ";
25         }
26         cout << endl;
27     }
28     cout << endl;
29     for (int i = 0; i < 10; ++i)
30     {
31         for (int j = 0; j <18; j++)
32         {
33             cout << zout[i *18 + j] << " ";
34         }
35         cout << endl;
36     }
37  return 0;
38 }