/*                   Kriging Interpolator                                            written by Chao-yi Lang                     July, 1995                     lang@cs.cornell.edu*/ #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 );     }