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 }
