kriging.c

Kriging 插值程序

/*                   Kriging Interpolator                       
                     written by Chao-yi Lang
                     July, 1995
                     [e-mail]

*/
 
#include <dx/dx.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

char krig_fn[]="_krig.dx";
char vari_fn[]="_vari.dx";
float *D,*weight;
float result[2];

float kriging_result(float *V, float i_f, float j_f, int dim, float *pos ,
        float *Z_s, int mode,int a, float c0, float c1);
int kriging(int *range, int mode, int item, float *Z_s, int *resol,
            float *pos, float c0, float c1,float a);
int *get_range(float *pos, int *range, int item);


Error m_Kriging(Object *in, Object *out)
{
    Array a1,a2,a3,a4;
    Object o = NULL;
    float x,*pos,*value,nugget,range_a,sill;
    int i, *resol_a, *range, range_default[4];
    char *method_a;
    char method_default[] = "exponential";

    
    int item, n, n1, rank, rank1, shape[3], shape1[3], method;
    Category category;
    Type type;

    /* copy the structure of in[0] */

    if (!in[0])  {
          DXSetError(ERROR_BAD_PARAMETER, "missing data parameter");
          goto error;
    }

    if ( DXGetObjectClass(in[0]) != CLASS_FIELD) {
          DXSetError(ERROR_BAD_PARAMETER, "data is not a field");
          goto error;
    }

    a1 = (Array) DXGetComponentValue((Field) in[0], "data");
    if (!a1) {
          DXSetError(ERROR_MISSING_DATA, "field has no data");
          return ERROR;
    }

    DXGetArrayInfo(a1, &n, &type, &category, &rank, shape);

    value = (float *) DXGetArrayData(a1);

    a1 = (Array) DXGetComponentValue((Field) in[0], "positions");
    if (!a1) {
          DXSetError(ERROR_MISSING_DATA, "field has no positions");
          return ERROR;
    }

    DXGetArrayInfo(a1, &n1, &type, &category, &rank1, shape1);

    if ((rank1==1) && (shape1[0]==2)) {
          pos = (float *) DXGetArrayData(a1);
    }
    else {
        DXSetError(ERROR_MISSING_DATA, "input field must be rank:1 shape:2");
        return ERROR;
    }

    if (!in[1]) {
	range = get_range( pos, range_default, n1*2 );
    }
    else {
        if ( DXGetObjectClass(in[1]) != CLASS_ARRAY) {
            DXSetError(ERROR_BAD_PARAMETER,"bad class");
            goto error;
        }

        a2 = (Array)DXCopy(in[1], COPY_STRUCTURE);
        if (!a2)
            goto error;

        if ( DXGetArrayClass(a2) != CLASS_ARRAY) {
            DXSetError(ERROR_BAD_PARAMETER,"bad array class");
            goto error;
        }
        DXGetArrayInfo(a2,&item,&type,&category,&rank,shape);

        if (item != 2) {
            DXSetError(ERROR_BAD_PARAMETER, "range error");
            goto error;
        }
        range = (int *) DXGetArrayData(a2);
    }

    if (!in[2]) {
        DXSetError(ERROR_BAD_PARAMETER, "missing data parameter");
        goto error;
    }
    if ( DXGetObjectClass(in[2]) != CLASS_ARRAY) {
        DXSetError(ERROR_BAD_PARAMETER,"resolution error");
        goto error;
    }

    a1 = (Array)DXCopy(in[2], COPY_STRUCTURE);

    if (!a1)
        goto error;

    DXGetArrayInfo(a1,&item,&type,&category,&rank,shape);

    resol_a = (int *) DXGetArrayData(a1);

    if (item != 1) {
        DXSetError(ERROR_BAD_PARAMETER, "resolution error");
        goto error;
    }


    if (!in[3]) 
        method_a = method_default;
    else if ( DXGetObjectClass(in[3]) == CLASS_STRING ) 
    	    DXExtractString(in[3], &method_a);
    	 else {
            DXSetError(ERROR_BAD_PARAMETER,"method should be string");
            goto error;
    	 }

    if ( strcasecmp( method_a, "spherical") == 0
       || strcasecmp( method_a, "1") == 0)
        method = 1;
    else if ( strcasecmp( method_a, "exponential") == 0
       || strcasecmp( method_a, "2") == 0)
        method = 2;
    else if ( strcasecmp( method_a, "gaussian") == 0
       || strcasecmp( method_a, "3") == 0)
        method = 3;
    else {
        DXSetError(ERROR_BAD_PARAMETER,"input method error");
        goto error;
    }
    
    if (!in[4])  {
        DXSetError(ERROR_BAD_PARAMETER, "missing data parameter");
        goto error;
    }
    else 
        DXExtractFloat(in[4], &range_a);
         
    if (!in[5])  {
        DXSetError(ERROR_BAD_PARAMETER, "missing data parameter");
        goto error;
    }
    else 
        DXExtractFloat(in[5], &nugget);
         
    if (!in[6]) { 
        DXSetError(ERROR_BAD_PARAMETER, "missing data parameter");
        goto error;
    }
    else 
        DXExtractFloat(in[6], &sill); 

    kriging(range, method, n1, value, resol_a, pos, nugget, sill-nugget,
            range_a);
    
    out[0]=DXImportDX(krig_fn,NULL,NULL,NULL,NULL);
    out[1]=DXImportDX(vari_fn,NULL,NULL,NULL,NULL);
    unlink(krig_fn);
    unlink(vari_fn);
    return OK;

    error:
    return ERROR;
}




int kriging(int *range, int mode, int item, float *Z_s, int *resol,
            float *pos, float c0, float c1, float a) {

    FILE *opt, *opt1;
    int dim,i,j,*ipiv,info;
    float xs,ys,zs,i_f,j_f,k,l,cnt_x,cnt_y;
    float begin_row,begin_col,end_row,end_col,*V,*Cd,*work,test_t;
    int resol_x, resol_y;


    /* initialize values */
    begin_row = (float) *(range);
    begin_col = (float) *(range+1);
    end_row = (float) *(range+2);
    end_col = (float) *(range+3);
    dim = item + 1;
    resol_x = *(resol);
    resol_y = *(resol+1);
    

    /* allocate V array */
    V = (float *) DXAllocate (sizeof (float) * dim * dim );
    D = (float *) DXAllocate (sizeof (float) * dim);
    weight = (float *) DXAllocate (sizeof (float) * dim);

    /* allocate Cd array */
    Cd = (float *) DXAllocate (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]= sqrt( test_t );
        }
    for ( i=0; i< dim-1 ;i++) {
        V[i*dim+dim-1]= 1;
        V[(dim-1)*(dim)+i]=1;
    }
    V[(dim-1)*(dim)+i] = 0;

    /* caculate the variogram of sample datas and put into  V array */
    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] = 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] = c0 + c1 *( 1 - 
                                  exp(-3*Cd[i*dim+j]/a) );
                         break;
                case 3 : /* Gauss mode */
                         V[i*dim+j] = V[j*dim+i] = c0 + c1 *( 1 -
                                  exp(-3*Cd[i*dim+j]*Cd[i*dim+j]/a/a));
                         break;
                default: fprintf(stderr, " mode error \n");
                         break;
            }
        }

    /* release Cd array */
    DXFree(Cd);
    ipiv = (int *) DXAllocate (sizeof(int)*dim);
    work = (float *) DXAllocate (sizeof(float)*dim);

    /* call inverse matrix function to inverse matrix C */
    sgetrf_(&dim,&dim,V,&dim,ipiv,&info);
    if ( info != 0) fprintf(stderr, "matrix error!!\n");
    sgetri_(&dim,V,&dim,ipiv,work,&dim,&info);
    if ( info != 0) fprintf(stderr, "matrix error!!\n");

    DXFree(ipiv);
    DXFree(work);

    /* create a temperate file to put result */
    opt = fopen(krig_fn,"w");
    if ( opt == NULL ) {
        fprintf(stderr,"error open output file\n");
        return(-1);
    }

    opt1 = fopen(vari_fn,"w");
    if ( opt1 == NULL ) {
        fprintf(stderr,"error open output file\n");
        return(-1);
    }

    /* for loop for each point of the estimated block */
    cnt_x = (end_row - begin_row)/ (float) resol_x;
    cnt_y = (end_col - begin_col)/ (float) resol_y;
    fprintf(opt,"object 1 class gridconnections counts %d %d\n",
            resol_x+1,resol_y+1);
    fprintf(opt,"attribute \"element type\" string \"quads\"\n");
    fprintf(opt,"attribute \"ref\" string \"positions\"\n");
    fprintf(opt,"#\n");
    fprintf(opt,"object 2 class array type float rank 1 shape 2 items %d data follows\n",(resol_x+1)*(resol_y+1));

    fprintf(opt1,"object 1 class gridconnections counts %d %d\n",
            resol_x+1,resol_y+1);
    fprintf(opt1,"attribute \"element type\" string \"quads\"\n");
    fprintf(opt1,"attribute \"ref\" string \"positions\"\n");
    fprintf(opt1,"#\n");
    fprintf(opt1,"object 2 class array type float rank 1 shape 2 items %d data follows\n",(resol_x+1)*(resol_y+1));

    for ( i = 0; i<= resol_x; i++) {
        i_f = cnt_x*i+begin_row; 
        for ( j = 0; j<= resol_y; j++) {
           j_f=cnt_y*j+begin_col;
          /* call kriging subroutine */
           fprintf(opt,"%f %f\n",i_f,j_f);
           fprintf(opt1,"%f %f\n",i_f,j_f);
	   l++;
        }
    }
    fprintf(opt,"#\n");
    fprintf(opt1,"#\n");
    fprintf(opt,"object 3 class array type float rank 0 items %d data follows\n", (resol_x+1)*(resol_y+1));
    fprintf(opt1,"object 3 class array type float rank 0 items %d data follows\n", (resol_x+1)*(resol_y+1));
    for ( i = 0; i<= resol_x; i++) {
        i_f = cnt_x*i+begin_row; 
        for ( j = 0; j<= resol_y; j++) {
           j_f=cnt_y*j+begin_col;
          /* call kriging subroutine */
           kriging_result(V,i_f,j_f,dim,pos,Z_s,mode,a,c0,c1);
           fprintf(opt,"%f\n",result[0]);
           fprintf(opt1,"%f\n",result[1]);
        }
    }
    fprintf(opt,"attribute \"dep\" string \"positions\"\n");
    fprintf(opt,"#\n");
    fprintf(opt,"object \"default\" class field\n");
    fprintf(opt,"component \"connections\" value 1\n");
    fprintf(opt,"component \"positions\" value 2\n");
    fprintf(opt,"component \"data\" value 3\n");
    fprintf(opt,"#\n");
    fprintf(opt,"end\n");
    fclose(opt);
    fprintf(opt1,"attribute \"dep\" string \"positions\"\n");
    fprintf(opt1,"#\n");
    fprintf(opt1,"object \"default\" class field\n");
    fprintf(opt1,"component \"connections\" value 1\n");
    fprintf(opt1,"component \"positions\" value 2\n");
    fprintf(opt1,"component \"data\" value 3\n");
    fprintf(opt1,"#\n");
    fprintf(opt1,"end\n");
    fclose(opt1);
    DXFree(V);
    DXFree(D);
    DXFree(weight);
    return OK;
error:
    return ERROR;
}


/*--------- kriging_result subroutine -----------*/
float kriging_result(float *V, float i_f, float j_f, int dim, float *pos,
        float *Z_s, int mode,int a, float c0, float c1) {

    int i,j;
    float h;

/* caculate the distance between estimated point and sample datas */
/* and caculate the variogram of estimated point and sample datas and
          put into D array */
    for (i=0; i < dim-1 ; i++) {
        h = sqrt( ( pos[i*2]-i_f )*( pos[i*2]-i_f ) +
                  ( pos[i*2+1]-j_f )*( pos[i*2+1]-j_f ) );
        switch( mode ) {
            case 1 : /* Spher mode */
                     if ( h < a )
                        D[i] = c0 + c1 * (1.5*h/a - 0.5*(h/a)*(h/a)*(h/a));
                     else
                        D[i] = c0 + c1;
                     break;
            case 2 : /* Expon mode */
                     D[i] = c0 + c1 * (1 - exp(-3*h/a));
                     break;
            case 3 : /* Gauss mode */
                     D[i] = c0 + c1 * (1 - exp(-3*h*h/a/a));
                     break;
            default: fprintf(stderr, " mode error \n");
                     break;
        }
    }
    D[i] = 1;

    /* caculate the weights */
    for ( i = 0 ;i <= dim ;i++) {
       weight[i] = 0;
       for ( j = 0 ; j <= dim ; j++)
          weight[i] += V[i*dim+j] * D[j];
    }

    /* caculate and return the estimated value */
    result[0] = 0;
    for ( i = 0 ;i < dim ;i++) {
        result[0] += weight[i] * Z_s[i];
    }
/* result[0] : kriging result */
/* result[1] : error variance result */
    result[0] = ( result[0] >= 0 ) ? result[0] : 0;

    result[1] = 0;
    for ( i = 0 ;i < dim-1 ;i++) {
        result[1] += weight[i] * D[i];
    }
    result[1] += weight[dim-1];
    result[1] = sqrt(result[1]);
    return (1);
}


/* get the maximum and minimum value of input range as default range */
int *get_range(float *pos, int *range, int item) {

    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=2; i<item; i+=2){
        if ( min_row > *(pos+i) ) min_row = *(pos+i);
        if ( min_col > *(pos+i+1) ) min_col = *(pos+i+1);
        if ( max_row < *(pos+i) ) max_row = *(pos+i);        
        if ( max_col < *(pos+i+1) ) max_col = *(pos+i+1);
    }
    
    *range = (int) min_row;
    *(range+1) = (int) min_col;
    *(range+2) = (int) max_row+1;      
    *(range+3) = (int) max_col+1;
    
    return ( range ); 
    
}