正文
C语言实现matlab的interp2()函数
提示：扫一扫查出行【扫一扫了解最新限行尾号】
复制提示
项目要用到matlab中的Vq = interp2(X,Y,V,Xq,Yq)函数，即把一个已知经纬度和对应值的矩阵，插值变换到一个给定经纬度网格中，也就是对给定网格填值，需要用到插值，这里使用双线性内插法。
*(这只是一个初步完成代码，仅供参考)
以下是对应C代码和测试程序：
//************************************
// 函数名称: inter_linear()
// 函数说明：计算两点之间某一给定点的值，（x0,y0）->(x1,y1),已知x(x0<x<x1),求y
// 返 回 值: double
// 参 数: double x0,y0,x1,y1/*in*/ 两个点的坐标（x0,y0）（x1,y1） // 作 者：WSS
// 作成日期：2019/09/05 //************************************
double inter_linear(double x0, double y0, double x1, double y1, double x)
{
double a0, a1, y;
a0 = (x - x1) / (x0 - x1);
a1 = (x - x0) / (x1 - x0);
y = a0*y0 + a1*y1;
return y;
}
//************************************
// 函数名称: interp2d()
// 函数说明：二维插值，同matlab的interp2()功能
// 返 回 值: double
// 参 数: x,y分别为长度为m和n的向量（一维数组），z为矩阵（对应的二维数组(m,n)）
// a,b分别为长度为asize和bsize的向量（一维数组），out_result为矩阵（对应的二维数组(asize,bsize)） // 作 者：WSS
// 作成日期：2019/09/05 //************************************
int interp2d(double *x, double *y, double *z, int m, int n, double *a, double *b, int asize,int bsize,double *out_result)
{
//a,b是待插值的矩阵行和列，遍历a，取出行，遍历b取出列，将插值计算后的值存储在out_result中
double pointa = ;
double pointb = ;
double w1, w2,w;
const int nu_val = -;
int tempi = ;
int tempj =;
for (int i = ; i < asize; ++i)
{
for (int j = ; j < bsize; ++j)
{
// 找到网格点位置
pointa = a[i];
pointb = b[j];
//确定pointa pointb所在的位置
for (int p = ; p < m-; ++p)
{
if (pointa < x[]||pointa>x[m-])//范围外的值不进行插值处理
{
tempi = nu_val;
break;
}
else if (pointa == x[m - ])
{
tempi = m - ;
break;
}
else if(pointa >= x[p] && pointa < x[p+])//x升序
{
tempi = p;
break;
}
}
for (int q = ; q < n - ; ++q)
{
if (pointb< y[] || pointb>y[n- ])//范围外的值不进行插值处理
{
tempj = nu_val;
break;
}
else if (pointb == y[n - ])
{
tempj = n - ;
break;
}
else if (pointb >= y[q] && pointb < y[q + ])//y升序
{
tempj = q;
break;
}
}
if (tempj == nu_val || tempi == nu_val)
{
out_result[i*bsize + j] = nu_val;
}
else
{
//如果四个点，有一个点为-9999，就不进行插值了
if(z[tempi*n + tempj]<-||z[tempi*n + tempj + ]<-||z[(tempi + )*n + tempj]<-||z[(tempi + )*n + tempj + ]<-)
{
out_result[i*bsize + j] = nu_val;
}
else
{
//cout << tempi << " " << tempj << endl;
/**************x方向进行插值*************/
if (x[tempi] == pointa)
{
//取网格节点值
w1 = z[tempi*n + tempj];
w2 = z[tempi*n + tempj + ];
//y方向进行插值
if (y[tempj] == pointb)
{
w = w1;
}
else
{
w = inter_linear(y[tempj], w1, y[tempj + ], w2, pointb);
}
}
else
{
//x方向进行插值
w1 = inter_linear(x[tempi], z[tempi*n + tempj], x[tempi + ], z[(tempi + )*n + tempj], pointa);
w2 = inter_linear(x[tempi], z[tempi*n + tempj + ], x[tempi + ], z[(tempi + )*n + tempj + ], pointa);
/*******y方向进行插值********/ if (y[tempj] == pointb) { //取网格节点值 w = w1; }
else
{
//进行插值(y) w = inter_linear(y[tempj], w1, y[tempj + ], w2, pointb);
}
}
out_result[i*bsize + j] = w;
} } }
}
return ;
}
测试函数：
int main()
{
//插值函数测试
int ret = ;
double a1[] = { , , , , };
double a2[] = {,,,,,};
double z[] ;
for (int i = ; i < ; ++i)
{
for (int j = ; j < ; ++j)
{
z[i*+j] = a1[i]+a2[j];
}
} double b1[] = { , , , , , , , , , };
double b2[] = { , , , , , , , , , , , , , , , , , };
double zout[];
ret = interp2d(a1, a2, z, , , b1, b2, , , zout);
for (int i = ; i < ; ++i)
{
for (int j = ; j < ; j++)
{
cout << z[i * + j] << " ";
}
cout << endl;
}
cout << endl;
for (int i = ; i < ; ++i)
{
for (int j = ; j <; j++)
{
cout << zout[i * + j] << " ";
}
cout << endl;
}
return 0;
}
