GRASS Programmer's Manual  6.4.2(2012)
interp2d.c
Go to the documentation of this file.
00001 
00002 /*-
00003  *
00004  * Original program and various modifications:
00005  * Lubos Mitas 
00006  *
00007  * GRASS4.1 version of the program and GRASS4.2 modifications:
00008  * H. Mitasova,
00009  * I. Kosinovsky, D. Gerdes
00010  * D. McCauley 
00011  *
00012  * Copyright 1993, 1995:
00013  * L. Mitas ,
00014  * H. Mitasova,
00015  * I. Kosinovsky,
00016  * D.Gerdes 
00017  * D. McCauley
00018  *
00019  * modified by McCauley in August 1995
00020  * modified by Mitasova in August 1995, Nov. 1996
00021  * bug fixes(mask) and modif. for variable smoothing Mitasova Jan 1997
00022  *
00023  */
00024 
00025 
00026 #include <stdio.h>
00027 #include <math.h>
00028 #include <unistd.h>
00029 #include <grass/gis.h>
00030 #include <grass/glocale.h>
00031 #include <grass/bitmap.h>
00032 
00033 #include <grass/interpf.h>
00034 
00035 
00036 #define CEULER .57721566
00037 
00038 
00039 int IL_grid_calc_2d(struct interp_params *params, struct quaddata *data,        /* given segment */
00040                     struct BM *bitmask, /* bitmask */
00041                     double zmin, double zmax,   /* min and max input z-values */
00042                     double *zminac, double *zmaxac,     /* min and max interp. z-values */
00043                     double *gmin, double *gmax, /* min and max inperp. slope val. */
00044                     double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */
00045                     double *ertot,      /* total interplating func. error */
00046                     double *b,  /* solutions of linear equations */
00047                     int offset1,        /* offset for temp file writing */
00048                     double dnorm)
00049 
00050 /*
00051  * Calculates grid for the given segment represented by data (contains
00052  * n_rows, n_cols, ew_res,ns_res, and all points inside + overlap) using
00053  * solutions of system of lin. equations and interpolating functions
00054  * interp() and interpder(). Also calls secpar() to compute slope, aspect
00055  * and curvatures if required.
00056  */
00057 {
00058 
00059     /*
00060      * C C       INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
00061      * c
00062      */
00063     double x_or = data->x_orig;
00064     double y_or = data->y_orig;
00065     int n_rows = data->n_rows;
00066     int n_cols = data->n_cols;
00067     int n_points = data->n_points;
00068     struct triple *points;
00069     static double *w2 = NULL;
00070     static double *w = NULL;
00071     int cond1, cond2;
00072     double r;
00073     double stepix, stepiy, xx, xg, yg, xx2;
00074     double rfsta2, /* cons, cons1, */ wm, dx, dy, dxx, dyy, dxy, h, bmgd1,
00075         bmgd2;
00076     double r2, gd1, gd2;        /* for interpder() */
00077     int n1, k, l, m;
00078     int ngstc, nszc, ngstr, nszr;
00079     double zz;
00080     int bmask = 1;
00081     static int first_time_z = 1;
00082     int offset, offset2;
00083     double fstar2 = params->fi * params->fi / 4.;
00084     double tfsta2, tfstad;
00085     double ns_res, ew_res;
00086     double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
00087     double xxr, yyr;
00088 
00089     if (params->theta) {
00090         teta = params->theta / 57.295779;       /* deg to rad */
00091         rsin = sin(teta);
00092         rcos = cos(teta);
00093     }
00094     if (params->scalex)
00095         scale = params->scalex;
00096 
00097     ns_res = (((struct quaddata *)(data))->ymax -
00098               ((struct quaddata *)(data))->y_orig) / data->n_rows;
00099     ew_res = (((struct quaddata *)(data))->xmax -
00100               ((struct quaddata *)(data))->x_orig) / data->n_cols;
00101 
00102     /*  tfsta2 = fstar2 * 2.; modified after removing normalization of z */
00103     tfsta2 = (fstar2 * 2.) / dnorm;
00104     tfstad = tfsta2 / dnorm;
00105     points = data->points;
00106 
00107     /*
00108      * normalization
00109      */
00110     stepix = ew_res / dnorm;
00111     stepiy = ns_res / dnorm;
00112 
00113     cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
00114              (params->adxy != NULL));
00115     cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);
00116 
00117     if (!w) {
00118         if (!(w = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
00119             G_warning(_("Out of memory"));
00120             return -1;
00121         }
00122     }
00123     if (!w2) {
00124         if (!(w2 = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
00125             G_warning(_("Out of memory"));
00126             return -1;
00127         }
00128     }
00129     n1 = n_points + 1;
00130     /*
00131      * C C         INTERPOLATION   *  MOST INNER LOOPS ! C
00132      */
00133     ngstc = (int)(x_or / ew_res + 0.5) + 1;
00134     nszc = ngstc + n_cols - 1;
00135     ngstr = (int)(y_or / ns_res + 0.5) + 1;
00136     nszr = ngstr + n_rows - 1;
00137 
00138 
00139     for (k = ngstr; k <= nszr; k++) {
00140         offset = offset1 * (k - 1);     /* rows offset */
00141         yg = (k - ngstr) * stepiy + stepiy / 2.;        /* fixed by J.H. in July 01 */
00142         for (m = 1; m <= n_points; m++) {
00143             wm = yg - points[m - 1].y;
00144             w[m] = wm;
00145             w2[m] = wm * wm;
00146         }
00147         for (l = ngstc; l <= nszc; l++) {
00148             if (bitmask != NULL)
00149                 /*      if(params->maskmap != NULL)  PK Apr 03 MASK support */
00150                 bmask = BM_get(bitmask, l - 1, k - 1);  /*fixed by helena jan 97 */
00151             /*    if(bmask==0 || bmask==-1) fprintf(stderr, "bmask=%d, at (%d,%d)\n", bmask, l, k); */
00152             xg = (l - ngstc) * stepix + stepix / 2.;    /*fixed by J.H. in July 01 */
00153             dx = 0.;
00154             dy = 0.;
00155             dxx = 0.;
00156             dyy = 0.;
00157             dxy = 0.;
00158             zz = 0.;
00159             if (bmask == 1) {   /* compute everything for area which is
00160                                  * not masked out */
00161                 h = b[0];
00162                 for (m = 1; m <= n_points; m++) {
00163                     xx = xg - points[m - 1].x;
00164                     if ((params->theta) && (params->scalex)) {
00165                         /* we run anisotropy */
00166                         xxr = xx * rcos + w[m] * rsin;
00167                         yyr = w[m] * rcos - xx * rsin;
00168                         xx2 = xxr * xxr;
00169                         w2[m] = yyr * yyr;
00170                         r2 = scale * xx2 + w2[m];
00171                         r = r2;
00172                         rfsta2 = scale * xx2 + w2[m];
00173                     }
00174                     else {
00175                         xx2 = xx * xx;
00176                         r2 = xx2 + w2[m];
00177                         r = r2;
00178                         rfsta2 = xx2 + w2[m];
00179                     }
00180 
00181                     h = h + b[m] * params->interp(r, params->fi);
00182                     if (cond1) {
00183                         if (!params->interpder(r, params->fi, &gd1, &gd2))
00184                             return -1;
00185                         bmgd1 = b[m] * gd1;
00186                         dx = dx + bmgd1 * xx;
00187                         dy = dy + bmgd1 * w[m];
00188                         if (cond2) {
00189                             bmgd2 = b[m] * gd2;
00190                             dxx = dxx + bmgd2 * xx2 + bmgd1;
00191                             dyy = dyy + bmgd2 * w2[m] + bmgd1;
00192                             dxy = dxy + bmgd2 * xx * w[m];
00193                         }
00194                     }
00195                 }
00196 
00197                 /*      zz = (h * dnorm) + zmin; replaced by helena jan. 97 due to removing norma
00198                    lization of z and zm in segmen2d.c */
00199                 zz = h + zmin;
00200                 if (first_time_z) {
00201                     first_time_z = 0;
00202                     *zmaxac = *zminac = zz;
00203                 }
00204                 *zmaxac = amax1(zz, *zmaxac);
00205                 *zminac = amin1(zz, *zminac);
00206                 if ((zz > zmax + 0.1 * (zmax - zmin))
00207                     || (zz < zmin - 0.1 * (zmax - zmin))) {
00208                     static int once = 0;
00209 
00210                     if (!once) {
00211                         once = 1;
00212                         G_warning(_("Overshoot - increase in tension suggested. "
00213                                     "Overshoot occures at (%d,%d) cell. "
00214                                     "Z-value %f, zmin %f, zmax %f."),
00215                                   l, k, zz, zmin, zmax);
00216                     }
00217                 }
00218                 
00219                 params->az[l] = (FCELL) zz;
00220                 
00221                 if (cond1) {
00222                     params->adx[l] = (FCELL) (-dx * tfsta2);
00223                     params->ady[l] = (FCELL) (-dy * tfsta2);
00224                     if (cond2) {
00225                         params->adxx[l] = (FCELL) (-dxx * tfstad);
00226                         params->adyy[l] = (FCELL) (-dyy * tfstad);
00227                         params->adxy[l] = (FCELL) (-dxy * tfstad);
00228                     }
00229                 }
00230 
00231             }
00232             else {
00233                 G_set_d_null_value(params->az + l, 1);
00234                 /*          fprintf (stderr, "zz=%f, az[l]=%f, c=%d\n", zz, params->az[l], l); */
00235 
00236                 if (cond1) {
00237                     G_set_d_null_value(params->adx + l, 1);
00238                     G_set_d_null_value(params->ady + l, 1);
00239                     if (cond2) {
00240                         G_set_d_null_value(params->adxx + l, 1);
00241                         G_set_d_null_value(params->adyy + l, 1);
00242                         G_set_d_null_value(params->adxy + l, 1);
00243                     }
00244                 }
00245             }
00246 
00247         }
00248         if (cond1 && (params->deriv != 1)) {
00249             if (params->secpar(params, ngstc, nszc, k, bitmask,
00250                                gmin, gmax, c1min, c1max, c2min, c2max, cond1,
00251                                cond2) < 0)
00252                 return -1;
00253         }
00254 
00255         offset2 = (offset + ngstc - 1) * sizeof(FCELL);
00256         if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
00257             return -1;
00258     }
00259     return 1;
00260 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines