Kriging法(克里金法)算法 C#实现?

这个是用以C++写 里有详细的注释 可以参考

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 );

}