GRASS Programmer's Manual  6.4.2(2012)
segmen2d.c
Go to the documentation of this file.
00001 /*
00002  **  Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993 
00003  **  Copyright  H. Mitasova, I. Kosinovsky, D.Gerdes  
00004  */
00005 
00006 #include <stdio.h>
00007 #include <stdlib.h>
00008 #include <math.h>
00009 #include <grass/gis.h>
00010 #include <grass/glocale.h>
00011 #include <grass/interpf.h>
00012 #include <grass/gmath.h>
00013 
00014 static double smallest_segment(struct multtree *, int);
00015 
00016 int IL_interp_segments_2d(struct interp_params *params, struct tree_info *info, /* info for the quad tree */
00017                           struct multtree *tree,        /* current leaf of the quad tree */
00018                           struct BM *bitmask,   /* bitmask */
00019                           double zmin, double zmax,     /* min and max input z-values */
00020                           double *zminac, double *zmaxac,       /* min and max interp. z-values */
00021                           double *gmin, double *gmax,   /* min and max inperp. slope val. */
00022                           double *c1min, double *c1max, double *c2min, double *c2max,   /* min and max interp. curv. val. */
00023                           double *ertot,        /* total interplating func. error */
00024                           int totsegm,  /* total number of segments */
00025                           int offset1,  /* offset for temp file writing */
00026                           double dnorm)
00027 /*
00028    Recursively processes each segment in a tree by
00029    a) finding points from neighbouring segments so that the total number of
00030    points is between KMIN and KMAX2 by calling tree function MT_get_region().
00031    b) creating and solving the system of linear equations using these points
00032    and interp() by calling matrix_create() and G_ludcmp().
00033    c) checking the interpolating function values at points by calling
00034    check_points().
00035    d) computing grid for this segment using points and interp() by calling
00036    grid_calc().
00037  */
00038 {
00039     double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2;
00040     int i, npt, nptprev, MAXENC;
00041     struct quaddata *data;
00042     static int cursegm = 0;
00043     static double *b = NULL;
00044     static int *indx = NULL;
00045     static double **matrix = NULL;
00046     double ew_res, ns_res;
00047     static int first_time = 1;
00048     static double smseg;
00049     int MINPTS;
00050     double pr;
00051     struct triple *point;
00052     struct triple skip_point;
00053     int m_skip, skip_index, j, k, segtest;
00054     double xx, yy, zz;
00055 
00056     /* find the size of the smallest segment once */
00057     if (first_time) {
00058         smseg = smallest_segment(info->root, 4);
00059         first_time = 0;
00060     }
00061     ns_res = (((struct quaddata *)(info->root->data))->ymax -
00062               ((struct quaddata *)(info->root->data))->y_orig) /
00063         params->nsizr;
00064     ew_res =
00065         (((struct quaddata *)(info->root->data))->xmax -
00066          ((struct quaddata *)(info->root->data))->x_orig) / params->nsizc;
00067 
00068     if (tree == NULL)
00069         return -1;
00070     if (tree->data == NULL)
00071         return -1;
00072     if (((struct quaddata *)(tree->data))->points == NULL) {
00073         for (i = 0; i < 4; i++) {
00074             IL_interp_segments_2d(params, info, tree->leafs[i],
00075                                   bitmask, zmin, zmax, zminac, zmaxac, gmin,
00076                                   gmax, c1min, c1max, c2min, c2max, ertot,
00077                                   totsegm, offset1, dnorm);
00078         }
00079         return 1;
00080     }
00081     else {
00082         distx = (((struct quaddata *)(tree->data))->n_cols * ew_res) * 0.1;
00083         disty = (((struct quaddata *)(tree->data))->n_rows * ns_res) * 0.1;
00084         distxp = 0;
00085         distyp = 0;
00086         xmn = ((struct quaddata *)(tree->data))->x_orig;
00087         xmx = ((struct quaddata *)(tree->data))->xmax;
00088         ymn = ((struct quaddata *)(tree->data))->y_orig;
00089         ymx = ((struct quaddata *)(tree->data))->ymax;
00090         i = 0;
00091         MAXENC = 0;
00092         /* data is a window with zero points; some fields don't make sence in this case
00093            so they are zero (like resolution,dimentions */
00094         /* CHANGE */
00095         /* Calcutaing kmin for surrent segment (depends on the size) */
00096 
00097 /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
00098         pr = pow(2., (xmx - xmn) / smseg - 1.);
00099         MINPTS =
00100             params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
00101         /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
00102 
00103         data =
00104             (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
00105                                              xmx + distx, ymx + disty, 0, 0,
00106                                              0, params->KMAX2);
00107         npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
00108 
00109         while ((npt < MINPTS) || (npt > params->KMAX2)) {
00110             if (i >= 70) {
00111                 G_warning(_("Taking too long to find points for interpolation - "
00112                             "please change the region to area where your points are. "
00113                             "Continuing calculations..."));
00114                 break;
00115             }
00116             i++;
00117             if (npt > params->KMAX2)
00118                 /* decrease window */
00119             {
00120                 MAXENC = 1;
00121                 nptprev = npt;
00122                 temp1 = distxp;
00123                 distxp = distx;
00124                 distx = distxp - fabs(distx - temp1) * 0.5;
00125                 temp2 = distyp;
00126                 distyp = disty;
00127                 disty = distyp - fabs(disty - temp2) * 0.5;
00128                 /* decrease by 50% of a previous change in window */
00129             }
00130             else {
00131                 nptprev = npt;
00132                 temp1 = distyp;
00133                 distyp = disty;
00134                 temp2 = distxp;
00135                 distxp = distx;
00136                 if (MAXENC) {
00137                     disty = fabs(disty - temp1) * 0.5 + distyp;
00138                     distx = fabs(distx - temp2) * 0.5 + distxp;
00139                 }
00140                 else {
00141                     distx += distx;
00142                     disty += disty;
00143                 }
00144                 /* decrease by 50% of extra distance */
00145             }
00146             data->x_orig = xmn - distx; /* update window */
00147             data->y_orig = ymn - disty;
00148             data->xmax = xmx + distx;
00149             data->ymax = ymx + disty;
00150             data->n_points = 0;
00151             npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
00152         }
00153         
00154         if (totsegm != 0) {
00155             G_percent(cursegm, totsegm, 1);
00156         }
00157         data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
00158         data->n_cols = ((struct quaddata *)(tree->data))->n_cols;
00159 
00160         /* for printing out overlapping segments */
00161         ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
00162         ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
00163         ((struct quaddata *)(tree->data))->xmax = xmx + distx;
00164         ((struct quaddata *)(tree->data))->ymax = ymx + disty;
00165 
00166         data->x_orig = xmn;
00167         data->y_orig = ymn;
00168         data->xmax = xmx;
00169         data->ymax = ymx;
00170 
00171         if (!matrix) {
00172             if (!
00173                 (matrix =
00174                  G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
00175                 G_warning(_("Out of memory"));
00176                 return -1;
00177             }
00178         }
00179         if (!indx) {
00180             if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
00181                 G_warning(_("Out of memory"));
00182                 return -1;
00183             }
00184         }
00185         if (!b) {
00186             if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
00187                 G_warning(_("Out of memory"));
00188                 return -1;
00189             }
00190         }
00191         /* allocate memory for CV points only if cv is performed */
00192         if (params->cv) {
00193             if (!
00194                 (point =
00195                  (struct triple *)G_malloc(sizeof(struct triple) *
00196                                            data->n_points))) {
00197                 G_warning(_("Out of memory"));
00198                 return -1;
00199             }
00200         }
00201 
00202         /*normalize the data so that the side of average segment is about 1m */
00203         /* put data_points into point only if CV is performed */
00204 
00205         for (i = 0; i < data->n_points; i++) {
00206             data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
00207             data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
00208             if (params->cv) {
00209                 point[i].x = data->points[i].x; /*cv stuff */
00210                 point[i].y = data->points[i].y; /*cv stuff */
00211                 point[i].z = data->points[i].z; /*cv stuff */
00212             }
00213 
00214             /* commented out by Helena january 1997 as this is not necessary
00215                although it may be useful to put normalization of z back? 
00216                data->points[i].z = data->points[i].z / dnorm;
00217                this made smoothing self-adjusting  based on dnorm
00218                if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
00219              */
00220         }
00221 
00222         /* cv stuff */
00223         if (params->cv)
00224             m_skip = data->n_points;
00225         else
00226             m_skip = 1;
00227 
00228         /* remove after cleanup - this is just for testing */
00229         skip_point.x = 0.;
00230         skip_point.y = 0.;
00231         skip_point.z = 0.;
00232 
00233         for (skip_index = 0; skip_index < m_skip; skip_index++) {
00234             if (params->cv) {
00235                 segtest = 0;
00236                 j = 0;
00237                 xx = point[skip_index].x * dnorm + data->x_orig +
00238                     params->x_orig;
00239                 yy = point[skip_index].y * dnorm + data->y_orig +
00240                     params->y_orig;
00241                 zz = point[skip_index].z;
00242                 if (xx >= data->x_orig + params->x_orig &&
00243                     xx <= data->xmax + params->x_orig &&
00244                     yy >= data->y_orig + params->y_orig &&
00245                     yy <= data->ymax + params->y_orig) {
00246                     segtest = 1;
00247                     skip_point.x = point[skip_index].x;
00248                     skip_point.y = point[skip_index].y;
00249                     skip_point.z = point[skip_index].z;
00250                     for (k = 0; k < m_skip; k++) {
00251                         if (k != skip_index && params->cv) {
00252                             data->points[j].x = point[k].x;
00253                             data->points[j].y = point[k].y;
00254                             data->points[j].z = point[k].z;
00255                             j++;
00256                         }
00257                     }
00258                 }               /* segment area test */
00259             }
00260             if (!params->cv) {
00261                 if (params->
00262                     matrix_create(params, data->points, data->n_points,
00263                                   matrix, indx) < 0)
00264                     return -1;
00265             }
00266             else if (segtest == 1) {
00267                 if (params->
00268                     matrix_create(params, data->points, data->n_points - 1,
00269                                   matrix, indx) < 0)
00270                     return -1;
00271             }
00272             if (!params->cv) {
00273                 for (i = 0; i < data->n_points; i++)
00274                     b[i + 1] = data->points[i].z;
00275                 b[0] = 0.;
00276                 G_lubksb(matrix, data->n_points + 1, indx, b);
00277         /* put here condition to skip ertot if not needed */
00278                 params->check_points(params, data, b, ertot, zmin, dnorm,
00279                                      skip_point);
00280             }
00281             else if (segtest == 1) {
00282                 for (i = 0; i < data->n_points - 1; i++)
00283                     b[i + 1] = data->points[i].z;
00284                 b[0] = 0.;
00285                 G_lubksb(matrix, data->n_points, indx, b);
00286                 params->check_points(params, data, b, ertot, zmin, dnorm,
00287                                      skip_point);
00288             }
00289         }                       /*end of cv loop */
00290 
00291         if (!params->cv)
00292             if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
00293                 (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
00294                 (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {
00295 
00296                 if (params->grid_calc(params, data, bitmask,
00297                                       zmin, zmax, zminac, zmaxac, gmin, gmax,
00298                                       c1min, c1max, c2min, c2max, ertot, b,
00299                                       offset1, dnorm) < 0)
00300                     return -1;
00301             }
00302 
00303         /* show after to catch 100% */
00304         cursegm++;
00305         if (totsegm < cursegm)
00306             G_debug(1, "%d %d", totsegm, cursegm);
00307         
00308         if (totsegm != 0) {
00309             G_percent(cursegm, totsegm, 1);
00310         }
00311         /* 
00312            G_free_matrix(matrix);
00313            G_free_ivector(indx);
00314            G_free_vector(b);
00315          */
00316         G_free(data->points);
00317         G_free(data);
00318     }
00319     return 1;
00320 }
00321 
00322 static double smallest_segment(struct multtree *tree, int n_leafs)
00323 {
00324     static int first_time = 1;
00325     int ii;
00326     static double minside;
00327     double side;
00328 
00329     if (tree == NULL)
00330         return 0;
00331     if (tree->data == NULL)
00332         return 0;
00333     if (tree->leafs != NULL) {
00334         for (ii = 0; ii < n_leafs; ii++) {
00335             side = smallest_segment(tree->leafs[ii], n_leafs);
00336             if (first_time) {
00337                 minside = side;
00338                 first_time = 0;
00339             }
00340             if (side < minside)
00341                 minside = side;
00342         }
00343     }
00344     else {
00345         side = ((struct quaddata *)(tree->data))->xmax -
00346             ((struct quaddata *)(tree->data))->x_orig;
00347         return side;
00348     }
00349     
00350     return minside;
00351 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines