📄 kriging.cpp
字号:
int WINAPI agaus(float *a,float *b,int n)
//解方程组ax=b,n为未知数个数,a为n*n矩阵,b为n*1矩阵
//结果放b里边
//采用了全主元消元法
{ int *js,l,k,i,j,is,p,q;
float d,t;
js=(int *)malloc(n*sizeof(int));
l=1;
for (k=0;k<=n-2;k++)
{ d=0.0;
for (i=k;i<=n-1;i++)
for (j=k;j<=n-1;j++)
{ t=(float)fabs(a[i*n+j]);
if (t>d) { d=t; js[k]=j; is=i;}
}
if (d+1.0==1.0) l=0;
else
{ if (js[k]!=k)
for (i=0;i<=n-1;i++)
{ p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if (is!=k)
{ for (j=k;j<=n-1;j++)
{ p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if (l==0)
{ free(js);
return(0);
}
d=a[k*n+k];
for (j=k+1;j<=n-1;j++)
{ p=k*n+j; a[p]=a[p]/d;}
b[k]=b[k]/d;
for (i=k+1;i<=n-1;i++)
{ for (j=k+1;j<=n-1;j++)
{ p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if (fabs(d)+1.0==1.0)
{ free(js);
return(0);
}
b[n-1]=b[n-1]/d;
for (i=n-2;i>=0;i--)
{ t=0.0;
for (j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for (k=n-1;k>=0;k--)
if (js[k]!=k)
{ t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
free(js);
return(1);
}
void kriging(int *range, int mode, int item, float *Z_s, int *resol,
float *pos, float c0, float c1, float a,float *result)
// *range所有特征点的Xmin,Ymin,Xmax,Ymax
// int mode, 计算半方差矩阵的模式 1 2 3三种情况
// int item,特征点数量
// float *Z_s,item个点的采样值
// int *resol,x方向和y方向的分辨率。即，resol_x为x方向的间隔数；
// resol_y为y方向的间隔数,最后总点数为：
//（resol_x+1）*（resol_y+1）
// float *pos:X1 Y1 X2 Y2 … Xitem Yitem 所有特征点的坐标
// float c0, float c1, float a 三个常量参数
// result 输出结果:（resol_x+1）*（resol_y+1）个点的插值结果
// x带表行位置，y代表列位置
{
int dim,i,j,k,l;
float i_f,j_f,cnt_x,cnt_y;
float begin_row,begin_col,end_row,end_col,*Cd,test_t;
int resol_x, resol_y;
float *D,*V,*temp;
FILE *fp; //定义文件用于存放结果
fp=fopen("test_out.txt","w");//打开文件
/* initialize values */
//得到范围坐标
begin_row = (float) *(range);
begin_col = (float) *(range+1);
end_row = (float) *(range+2);
end_col = (float) *(range+3);
//为特征点数+1
dim = item + 1;
//得到分辩率
resol_x = *(resol);
resol_y = *(resol+1);
/* allocate V D array */
//V为(n+1)*(n+1)矩阵，V*W=D
//D为(n+1)*1矩阵，V*W=D
//temp临时存放V，大小与V同
//W放在D中，故不为W分配空间
V = (float *)GlobalAllocPtr(GMEM_MOVEABLE,sizeof (float) * dim * dim);
D = (float *)GlobalAllocPtr(GMEM_MOVEABLE,sizeof (float) * dim );
temp = (float *)GlobalAllocPtr(GMEM_MOVEABLE,sizeof (float) * dim * dim);
/* allocate Cd array */
// Cd为(n+1)*(n+1)矩阵，用于存储距离
// n 即item
// D0 0 D0 1 … D0 item-1 1
// D1 0 D1 1 … D1 item-1 1
// Ditem-1 0 Ditem-1 1 … Ditem-1 item-1 1
// 1 1 1 0
Cd = (float *) GlobalAllocPtr(GMEM_MOVEABLE,sizeof (float) * dim * dim);
//计算上面的距离，用的欧氏方法
/* caculate the distance between sample datas put into Cd array*/
for ( i=0; i< dim-1 ;i++)
for (j=i; j < dim-1 ; j++)
{
test_t =( pos[i*2]-pos[j*2] )*( pos[i*2]-pos[j*2])+
( pos[i*2+1]-pos[j*2+1] )*( pos[i*2+1]-pos[j*2+1] );
Cd[i*dim+j]=(float) sqrt( test_t );
}
//补1
for ( i=0; i< dim-1 ;i++)
{
V[i*dim+dim-1]= 1;
V[(dim-1)*(dim)+i]=1;
}
//添0
V[(dim-1)*(dim)+i] = 0;
/* caculate the variogram of sample datas and put into V array */
//依据Dij计算采样点的半方差矩阵并放入V
for ( i=0; i< dim-1 ;i++)
for (j=i; j < dim-1; j++)
{//由于对称，只计算一半
switch( mode )
{//模式不同，公式不同
case 1 : /* Spher mode */
if ( Cd[i*dim+j] < a )
V[i*dim+j] = V[j*dim+i] = (float)(c0 + c1*(
1.5*Cd[i*dim+j]/a - 0.5*(Cd[i*dim+j]/a)*
(Cd[i*dim+j]/a)*(Cd[i*dim+j]/a)));
else
V[i*dim+j] = V[j*dim+i] = c0 + c1;
break;
case 2 : /* Expon mode */
V[i*dim+j] = V[j*dim+i] =(float)( c0 + c1 *( 1 -
exp(-3*Cd[i*dim+j]/a) ));
break;
case 3 : /* Gauss mode */
V[i*dim+j] = V[j*dim+i] = (float)(c0 + c1 *( 1 -
exp(-3*Cd[i*dim+j]*Cd[i*dim+j]/a/a)));
break;
default: V[i*dim+j] = V[j*dim+i] =(float)(1*Cd[i*dim+j]);
break;
}
}
/* release Cd array */
GlobalFreePtr(Cd);
//Cd完成任务，释放
//先暂存V,因GAUSS算法要改V
for(k=0;k<dim*dim;k++)
temp[k]=V[k];
//计算所有范围内，以分辨率为间隔的所有点
//例：range:
//2,4, 10,20
//resol:4,4 两方向均取4格，
//则待算行列方向步长为：cnt_x=（10-2）/4=2 cnt_y=(20-4)/4=4
//待算点为：
//2 4 2 8 2 12 2 16 2 20
//4 4 4 8 4 12 4 16 4 20
//6 4 6 8 6 12 6 16 6 20
//8 4 8 8 8 12 8 16 8 20
//10 4 10 8 10 12 10 16 10 20
/* for loop for each point of the estimated block */
cnt_x = (end_row - begin_row)/ (float) resol_x;//x方向步长
cnt_y = (end_col - begin_col)/ (float) resol_y;// y方向步长
l=0;
for ( i = 0; i<= resol_x; i++)
{
i_f = cnt_x*i+begin_row; //2 4 6 8 10
for ( j = 0; j<= resol_y; j++)
{
j_f=cnt_y*j+begin_col;//4 8 12 16 20
// 得到待计算点坐标：行：i_f
// 列：j_f
// 针对该点计算D的item个值,并补一个1
for(k=0;k<dim-1;k++)
{
test_t =( i_f-pos[2*k] ) *( i_f-pos[2*k] )+( j_f-pos[2*k+1])*( j_f-pos[2*k+1] );
test_t = (float)(sqrt( test_t ));
switch( mode )
{//模式不同，公式不同
case 1 : /* Spher mode */
if ( test_t < a )
D[k] = (float)(c0 + c1*(
1.5*test_t/a - 0.5*(test_t/a)*
(test_t/a)*(test_t/a)));
else
D[k] = c0 + c1;
break;
case 2 : /* Expon mode */
D[k] = (float)(c0 + c1 *( 1 -
exp(-3*test_t/a)));
break;
case 3 : /* Gauss mode */
D[k] = (float)(c0 + c1 *( 1 -
exp(-3*test_t*test_t/a/a)));
break;
default: D[k]=(float)(1*test_t);
break;
}
}
D[dim-1]=1;//并补一个1
//现在V有了,D有了,可解方程了
//先取回V
for(k=0;k<dim*dim;k++)
V[k]=temp[k];
agaus(V,D,dim);
//解出的结果:W1,W2,...,Witem,λ在D中，替换D原来的值
//下面计算该点的预测值,依据公式:
//W1Y1+W2Y2+...+WitemYitem有:
test_t=0;//用于累加
for(k=0;k<dim-1;k++)
test_t+=D[k]*Z_s[k];
*(result+l++)=test_t;
//得到该点的最终结果
//////////////////////////
//(int)(i_f+0.5) 坐标 x
//(int)(j_f+0.5) 坐标 y
//(int)(test_t+0.5) 值 z
//写入文件
fprintf(fp,"%d,%d,%d\n",(int)(i_f+0.5),(int)(j_f+0.5),(int)(test_t+0.5));
//////////////////////////
}
}
//关闭文件
fclose(fp);
GlobalFreePtr(V);
GlobalFreePtr(D);
GlobalFreePtr(temp);
}
int *get_range(float *pos, int *range, int item)
//item 点数
//所有点的坐标：行、列、行、列...
//range:返回左上角行，左上角列，右下角行，右下角列
{
float min_row, min_col, max_row, max_col;
int i;
min_row = *pos;
min_col = *(pos+1);
max_row = *pos;
max_col = *(pos+1);
for ( i=1; i<item; i++){
if ( min_row > *(pos+2*i) ) min_row = *(pos+2*i);
if ( min_col > *(pos+2*i+1) ) min_col = *(pos+2*i+1);
if ( max_row < *(pos+2*i) ) max_row = *(pos+2*i);
if ( max_col < *(pos+2*i+1) ) max_col = *(pos+2*i+1);
}
*range = (int) min_row;
*(range+1) = (int) min_col;
*(range+2) = (int) max_row;
*(range+3) = (int) max_col;
return ( range );
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -