江苏龙网

  1. 首页 > 生活 > >

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

http://ilnumerics.net/kriging-interpolation-net.html不谢。
■网友
有R的，要吗？
■网友
这个是用以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\u0026lt;=n-2;k++) { d=0.0; for (i=k;i\u0026lt;=n-1;i++) for (j=k;j\u0026lt;=n-1;j++) { t=(float)fabs(a); if (t\u0026gt;d) { d=t; js=j; is=i;} } if (d+1.0==1.0) l=0; else { if (js!=k) for (i=0;i\u0026lt;=n-1;i++) { p=i*n+k; q=i*n+js; t=a; a=a; a=t; } if (is!=k) { for (j=k;j\u0026lt;=n-1;j++) { p=k*n+j; q=is*n+j; t=a; a=a; a=t; } t=b; b=b; b=t; } } if (l==0) { free(js); return(0); } d=a; for (j=k+1;j\u0026lt;=n-1;j++) { p=k*n+j; a=a/d;} b=b/d; for (i=k+1;i\u0026lt;=n-1;i++) { for (j=k+1;j\u0026lt;=n-1;j++) { p=i*n+j; a=a-a*a; } b=b-a*b; } } d=a; if (fabs(d)+1.0==1.0) { free(js); return(0); } b=b/d; for (i=n-2;i\u0026gt;=0;i--) { t=0.0; for (j=i+1;j\u0026lt;=n-1;j++) t=t+a*b; b=b-t; } js=n-1; for (k=n-1;k\u0026gt;=0;k--) if (js!=k) { t=b; b=b]; b]=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\u0026lt; dim-1 ;i++) for (j=i; j \u0026lt; dim-1 ; j++) { test_t =( pos-pos )*( pos-pos)+ ( pos-pos )*( pos-pos ); Cd=(float) sqrt( test_t ); } //补1 for ( i=0; i\u0026lt; dim-1 ;i++) { V= 1; V=1; } //添0 V = 0; /* caculate the variogram of sample datas and put into V array */ //依据Dij计算采样点的半方差矩阵并放入V for ( i=0; i\u0026lt; dim-1 ;i++) for (j=i; j \u0026lt; dim-1; j++) {//由于对称，只计算一半 switch( mode ) {//模式不同，公式不同 case 1 : /* Spher mode */ if ( Cd \u0026lt; a ) V = V = (float)(c0 + c1*( 1.5*Cd/a - 0.5*(Cd/a)* (Cd/a)*(Cd/a))); else V = V = c0 + c1; break; case 2 : /* Expon mode */ V = V =(float)( c0 + c1 *( 1 - exp(-3*Cd/a) )); break; case 3 : /* Gauss mode */ V = V = (float)(c0 + c1 *( 1 - exp(-3*Cd*Cd/a/a))); break; default: V = V =(float)(1*Cd); break; } } /* release Cd array */ GlobalFreePtr(Cd); //Cd完成任务，释放 //先暂存V,因GAUSS算法要改V for(k=0;k\u0026lt;dim*dim;k++) temp=V; //计算所有范围内，以分辨率为间隔的所有点 //例：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\u0026lt;= resol_x; i++) { i_f = cnt_x*i+begin_row; //2 4 6 8 10 for ( j = 0; j\u0026lt;= 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\u0026lt;dim-1;k++) { test_t =( i_f-pos ) *( i_f-pos )+( j_f-pos)*( j_f-pos ); test_t = (float)(sqrt( test_t )); switch( mode ) {//模式不同，公式不同 case 1 : /* Spher mode */ if ( test_t \u0026lt; a ) D = (float)(c0 + c1*( 1.5*test_t/a - 0.5*(test_t/a)* (test_t/a)*(test_t/a))); else D = c0 + c1; break; case 2 : /* Expon mode */ D = (float)(c0 + c1 *( 1 - exp(-3*test_t/a))); break; case 3 : /* Gauss mode */ D = (float)(c0 + c1 *( 1 - exp(-3*test_t*test_t/a/a))); break; default: D=(float)(1*test_t); break; } } D=1;//并补一个1 //现在V有了,D有了,可解方程了 //先取回V for(k=0;k\u0026lt;dim*dim;k++) V=temp; agaus(V,D,dim); //解出的结果:W1,W2,...,Witem,λ在D中，替换D原来的值 //下面计算该点的预测值,依据公式: //W1Y1+W2Y2+...+WitemYitem有: test_t=0;//用于累加 for(k=0;k\u0026lt;dim-1;k++) test_t+=D*Z_s; *(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\",(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\u0026lt;item; i++){ if ( min_row \u0026gt; *(pos+2*i) ) min_row = *(pos+2*i); if ( min_col \u0026gt; *(pos+2*i+1) ) min_col = *(pos+2*i+1); if ( max_row \u0026lt; *(pos+2*i) ) max_row = *(pos+2*i); if ( max_col \u0026lt; *(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 ); }


推荐阅读


上一篇：旷视科技手机人脸解锁方案中的活体检测是咋做到有效防御人像翻拍攻击的?

下一篇：软件与计算机方面科班出身与非科班出身到底有没有区别