GRASS Programmer's Manual  6.4.2(2012)
point2d.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  *
00022  */
00023 
00024 
00025 #include <stdio.h>
00026 #include <math.h>
00027 #include <unistd.h>
00028 #include <grass/gis.h>
00029 #include <grass/site.h>
00030 #include <grass/Vect.h>
00031 #include <grass/dbmi.h>
00032 
00033 #define POINT2D_C
00034 #include <grass/interpf.h>
00035 
00036 int IL_check_at_points_2d(struct interp_params *params, struct quaddata *data,  /* current region */
00037                           double *b,    /* solution of linear equations */
00038                           double *ertot,        /* total error */
00039                           double zmin,  /* min z-value */
00040                           double dnorm, struct triple skip_point)
00041 
00042 /*
00043  * Checks if interpolating function interp() evaluates correct z-values at
00044  * given points. If smoothing is used calculate the maximum error caused
00045  * by smoothing.
00046  */
00047 {
00048     int n_points = data->n_points;      /* number of points */
00049     struct triple *points = data->points;       /* points for interpolation */
00050     double east = data->xmax;
00051     double west = data->x_orig;
00052     double north = data->ymax;
00053     double south = data->y_orig;
00054     double rfsta2, errmax, h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
00055     double skip_err;
00056     int n1, mm, m, cat;
00057     double fstar2;
00058     int inside;
00059 
00060     /*  Site *site; */
00061     char buf[1024];
00062 
00063 
00064     /*  if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
00065        G_fatal_error ("Memory error for site struct"); */
00066 
00067     fstar2 = params->fi * params->fi / 4.;
00068     errmax = .0;
00069     n1 = n_points + 1;
00070     for (mm = 1; mm <= n_points; mm++) {
00071         h = b[0];
00072         for (m = 1; m <= n_points; m++) {
00073             xx = points[mm - 1].x - points[m - 1].x;
00074             yy = points[mm - 1].y - points[m - 1].y;
00075             r2 = yy * yy + xx * xx;
00076             if (r2 != 0.) {
00077                 rfsta2 = fstar2 * r2;
00078                 r = r2;
00079                 h = h + b[m] * params->interp(r, params->fi);
00080             }
00081         }
00082         /* modified by helena january 1997 - normalization of z was
00083            removed from segm2d.c and interp2d.c
00084            hz = (h * dnorm) + zmin;
00085            zz = (points[mm - 1].z * dnorm) + zmin;
00086          */
00087         hz = h + zmin;
00088         zz = points[mm - 1].z + zmin;
00089         err = hz - zz;
00090         xmm = points[mm - 1].x * dnorm + params->x_orig + west;
00091         ymm = points[mm - 1].y * dnorm + params->y_orig + south;
00092         if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
00093             && (ymm >= south + params->y_orig) &&
00094             (ymm <= north + params->y_orig))
00095             inside = 1;
00096         else
00097             inside = 0;
00098 
00099         if (params->fddevi != NULL) {
00100 
00101             if (inside) {       /* if the point is inside the region */
00102                 Vect_reset_line(Pnts);
00103                 Vect_reset_cats(Cats2);
00104 
00105                 Vect_append_point(Pnts, xmm, ymm, zz);
00106                 cat = count;
00107                 Vect_cat_set(Cats2, 1, cat);
00108                 Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
00109 
00110                 db_zero_string(&sql2);
00111                 sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
00112                 db_append_string(&sql2, buf);
00113 
00114                 sprintf(buf, ", %f", err);
00115                 db_append_string(&sql2, buf);
00116                 db_append_string(&sql2, ")");
00117                 G_debug(3, db_get_string(&sql2));
00118 
00119                 if (db_execute_immediate(driver2, &sql2) != DB_OK) {
00120                     db_close_database(driver2);
00121                     db_shutdown_driver(driver2);
00122                     G_fatal_error("Cannot insert new row: %s",
00123                                   db_get_string(&sql2));
00124                 }
00125                 count++;
00126 
00127             }
00128         }
00129         (*ertot) += err * err;
00130     }
00131 
00132     /* cv stuff */
00133     if (params->cv) {
00134         h = b[0];
00135         for (m = 1; m <= n_points - 1; m++) {
00136             xx = points[m - 1].x - skip_point.x;
00137             yy = points[m - 1].y - skip_point.y;
00138             r2 = yy * yy + xx * xx;
00139             if (r2 != 0.) {
00140                 rfsta2 = fstar2 * r2;
00141                 r = r2;
00142                 h = h + b[m] * params->interp(r, params->fi);
00143             }
00144         }
00145         hz = h + zmin;
00146         zz = skip_point.z + zmin;
00147         skip_err = hz - zz;
00148         xmm = skip_point.x * dnorm + params->x_orig + west;
00149         ymm = skip_point.y * dnorm + params->y_orig + south;
00150 
00151         if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
00152             && (ymm >= south + params->y_orig) &&
00153             (ymm <= north + params->y_orig))
00154             inside = 1;
00155         else
00156             inside = 0;
00157 
00158         if (inside) {           /* if the point is inside the region */
00159             Vect_reset_line(Pnts);
00160             Vect_reset_cats(Cats2);
00161 
00162             Vect_append_point(Pnts, xmm, ymm, zz);
00163             cat = count;
00164             Vect_cat_set(Cats2, 1, cat);
00165             Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
00166 
00167             db_zero_string(&sql2);
00168             sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
00169             db_append_string(&sql2, buf);
00170 
00171             sprintf(buf, ", %f", skip_err);
00172             db_append_string(&sql2, buf);
00173             db_append_string(&sql2, ")");
00174             G_debug(3, db_get_string(&sql2));
00175 
00176             if (db_execute_immediate(driver2, &sql2) != DB_OK) {
00177                 db_close_database(driver2);
00178                 db_shutdown_driver(driver2);
00179                 G_fatal_error("Cannot insert new row: %s",
00180                               db_get_string(&sql2));
00181             }
00182             count++;
00183         }
00184     }                           /* cv */
00185 
00186 
00187     return 1;
00188 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines