GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
segmen2d.c
Go to the documentation of this file.
1/*!
2 * \file segmen2d.c
3 *
4 * \author H. Mitasova, I. Kosinovsky, D. Gerdes
5 *
6 * \copyright
7 * (C) 1993 by Helena Mitasova and the GRASS Development Team
8 *
9 * \copyright
10 * This program is free software under the
11 * GNU General Public License (>=v2).
12 * Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 */
16
17
18#include <stdio.h>
19#include <stdlib.h>
20#include <math.h>
21#include <grass/gis.h>
22#include <grass/glocale.h>
23#include <grass/interpf.h>
24#include <grass/gmath.h>
25
26
27/*!
28 * Interpolate recursively a tree of segments
29 *
30 * Recursively processes each segment in a tree by:
31 * - finding points from neighbouring segments so that the total number of
32 * points is between KMIN and KMAX2 by calling tree function MT_get_region().
33 * - creating and solving the system of linear equations using these points
34 * and interp() by calling matrix_create() and G_ludcmp().
35 * - checking the interpolating function values at points by calling
36 * check_points().
37 * - computing grid for this segment using points and interp() by calling
38 * grid_calc().
39 *
40 * \todo
41 * Isn't this in fact the updated version of the function (IL_interp_segments_new_2d)?
42 * The function IL_interp_segments_new_2d has the following, better behavior:
43 * The difference between this function and IL_interp_segments_2d() is making
44 * sure that additional points are taken from all directions, i.e. it finds
45 * equal number of points from neighboring segments in each of 8 neighborhoods.
46 */
48 struct tree_info *info, /*!< info for the quad tree */
49 struct multtree *tree, /*!< current leaf of the quad tree */
50 struct BM *bitmask, /*!< bitmask */
51 double zmin, double zmax, /*!< min and max input z-values */
52 double *zminac, double *zmaxac, /*!< min and max interp. z-values */
53 double *gmin, double *gmax, /*!< min and max inperp. slope val. */
54 double *c1min, double *c1max, /*!< min and max interp. curv. val. */
55 double *c2min, double *c2max, /*!< min and max interp. curv. val. */
56 double *ertot, /*!< total interplating func. error */
57 int totsegm, /*!< total number of segments */
58 off_t offset1, /*!< offset for temp file writing */
59 double dnorm)
60{
61 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2;
62 int i, npt, nptprev, MAXENC;
63 struct quaddata *data;
64 static int cursegm = 0;
65 static double *b = NULL;
66 static int *indx = NULL;
67 static double **matrix = NULL;
68 double ew_res, ns_res;
69 static int first_time = 1;
70 static double smseg;
71 int MINPTS;
72 double pr;
73 struct triple *point;
74 struct triple skip_point;
75 int m_skip, skip_index, j, k, segtest;
76 double xx, yy, zz;
77
78 /* find the size of the smallest segment once */
79 if (first_time) {
80 smseg = smallest_segment(info->root, 4);
81 first_time = 0;
82 }
83 ns_res = (((struct quaddata *)(info->root->data))->ymax -
84 ((struct quaddata *)(info->root->data))->y_orig) /
85 params->nsizr;
86 ew_res =
87 (((struct quaddata *)(info->root->data))->xmax -
88 ((struct quaddata *)(info->root->data))->x_orig) / params->nsizc;
89
90 if (tree == NULL)
91 return -1;
92 if (tree->data == NULL)
93 return -1;
94 if (((struct quaddata *)(tree->data))->points == NULL) {
95 for (i = 0; i < 4; i++) {
96 IL_interp_segments_2d(params, info, tree->leafs[i],
97 bitmask, zmin, zmax, zminac, zmaxac, gmin,
98 gmax, c1min, c1max, c2min, c2max, ertot,
99 totsegm, offset1, dnorm);
100 }
101 return 1;
102 }
103 else {
104 distx = (((struct quaddata *)(tree->data))->n_cols * ew_res) * 0.1;
105 disty = (((struct quaddata *)(tree->data))->n_rows * ns_res) * 0.1;
106 distxp = 0;
107 distyp = 0;
108 xmn = ((struct quaddata *)(tree->data))->x_orig;
109 xmx = ((struct quaddata *)(tree->data))->xmax;
110 ymn = ((struct quaddata *)(tree->data))->y_orig;
111 ymx = ((struct quaddata *)(tree->data))->ymax;
112 i = 0;
113 MAXENC = 0;
114 /* data is a window with zero points; some fields don't make sense in this case
115 so they are zero (like resolution,dimensions */
116 /* CHANGE */
117 /* Calcutaing kmin for surrent segment (depends on the size) */
118
119/*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
120 pr = pow(2., (xmx - xmn) / smseg - 1.);
121 MINPTS =
122 params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
123 /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
124
125 data =
126 (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
127 xmx + distx, ymx + disty, 0, 0,
128 0, params->KMAX2);
129 npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
130
131 while ((npt < MINPTS) || (npt > params->KMAX2)) {
132 if (i >= 70) {
133 G_warning(_("Taking too long to find points for interpolation - "
134 "please change the region to area where your points are. "
135 "Continuing calculations..."));
136 break;
137 }
138 i++;
139 if (npt > params->KMAX2)
140 /* decrease window */
141 {
142 MAXENC = 1;
143 nptprev = npt;
144 temp1 = distxp;
145 distxp = distx;
146 distx = distxp - fabs(distx - temp1) * 0.5;
147 temp2 = distyp;
148 distyp = disty;
149 disty = distyp - fabs(disty - temp2) * 0.5;
150 /* decrease by 50% of a previous change in window */
151 }
152 else {
153 nptprev = npt;
154 temp1 = distyp;
155 distyp = disty;
156 temp2 = distxp;
157 distxp = distx;
158 if (MAXENC) {
159 disty = fabs(disty - temp1) * 0.5 + distyp;
160 distx = fabs(distx - temp2) * 0.5 + distxp;
161 }
162 else {
163 distx += distx;
164 disty += disty;
165 }
166 /* decrease by 50% of extra distance */
167 }
168 data->x_orig = xmn - distx; /* update window */
169 data->y_orig = ymn - disty;
170 data->xmax = xmx + distx;
171 data->ymax = ymx + disty;
172 data->n_points = 0;
173 npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
174 }
175
176 if (totsegm != 0) {
177 G_percent(cursegm, totsegm, 1);
178 }
179 data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
180 data->n_cols = ((struct quaddata *)(tree->data))->n_cols;
181
182 /* for printing out overlapping segments */
183 ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
184 ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
185 ((struct quaddata *)(tree->data))->xmax = xmx + distx;
186 ((struct quaddata *)(tree->data))->ymax = ymx + disty;
187
188 data->x_orig = xmn;
189 data->y_orig = ymn;
190 data->xmax = xmx;
191 data->ymax = ymx;
192
193 if (!matrix) {
194 if (!
195 (matrix =
196 G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
197 G_warning(_("Out of memory"));
198 return -1;
199 }
200 }
201 if (!indx) {
202 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
203 G_warning(_("Out of memory"));
204 return -1;
205 }
206 }
207 if (!b) {
208 if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
209 G_warning(_("Out of memory"));
210 return -1;
211 }
212 }
213 /* allocate memory for CV points only if cv is performed */
214 if (params->cv) {
215 if (!
216 (point =
217 (struct triple *)G_malloc(sizeof(struct triple) *
218 data->n_points))) {
219 G_warning(_("Out of memory"));
220 return -1;
221 }
222 }
223
224 /*normalize the data so that the side of average segment is about 1m */
225 /* put data_points into point only if CV is performed */
226
227 for (i = 0; i < data->n_points; i++) {
228 data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
229 data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
230 if (params->cv) {
231 point[i].x = data->points[i].x; /*cv stuff */
232 point[i].y = data->points[i].y; /*cv stuff */
233 point[i].z = data->points[i].z; /*cv stuff */
234 }
235
236 /* commented out by Helena january 1997 as this is not necessary
237 although it may be useful to put normalization of z back?
238 data->points[i].z = data->points[i].z / dnorm;
239 this made smoothing self-adjusting based on dnorm
240 if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
241 */
242 }
243
244 /* cv stuff */
245 if (params->cv)
246 m_skip = data->n_points;
247 else
248 m_skip = 1;
249
250 /* remove after cleanup - this is just for testing */
251 skip_point.x = 0.;
252 skip_point.y = 0.;
253 skip_point.z = 0.;
254
255
256 /*** TODO: parallelize this loop instead of the LU solver! ***/
257 for (skip_index = 0; skip_index < m_skip; skip_index++) {
258 if (params->cv) {
259 segtest = 0;
260 j = 0;
261 xx = point[skip_index].x * dnorm + data->x_orig +
262 params->x_orig;
263 yy = point[skip_index].y * dnorm + data->y_orig +
264 params->y_orig;
265 zz = point[skip_index].z;
266 if (xx >= data->x_orig + params->x_orig &&
267 xx <= data->xmax + params->x_orig &&
268 yy >= data->y_orig + params->y_orig &&
269 yy <= data->ymax + params->y_orig) {
270 segtest = 1;
271 skip_point.x = point[skip_index].x;
272 skip_point.y = point[skip_index].y;
273 skip_point.z = point[skip_index].z;
274 for (k = 0; k < m_skip; k++) {
275 if (k != skip_index && params->cv) {
276 data->points[j].x = point[k].x;
277 data->points[j].y = point[k].y;
278 data->points[j].z = point[k].z;
279 j++;
280 }
281 }
282 } /* segment area test */
283 }
284 if (!params->cv) {
285 if (params->
286 matrix_create(params, data->points, data->n_points,
287 matrix, indx) < 0)
288 return -1;
289 }
290 else if (segtest == 1) {
291 if (params->
292 matrix_create(params, data->points, data->n_points - 1,
293 matrix, indx) < 0)
294 return -1;
295 }
296 if (!params->cv) {
297 for (i = 0; i < data->n_points; i++)
298 b[i + 1] = data->points[i].z;
299 b[0] = 0.;
300 G_lubksb(matrix, data->n_points + 1, indx, b);
301 /* put here condition to skip error if not needed */
302 params->check_points(params, data, b, ertot, zmin, dnorm,
303 skip_point);
304 }
305 else if (segtest == 1) {
306 for (i = 0; i < data->n_points - 1; i++)
307 b[i + 1] = data->points[i].z;
308 b[0] = 0.;
309 G_lubksb(matrix, data->n_points, indx, b);
310 params->check_points(params, data, b, ertot, zmin, dnorm,
311 skip_point);
312 }
313 } /*end of cv loop */
314
315 if (!params->cv)
316 if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
317 (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
318 (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {
319
320 if (params->grid_calc(params, data, bitmask,
321 zmin, zmax, zminac, zmaxac, gmin, gmax,
322 c1min, c1max, c2min, c2max, ertot, b,
323 offset1, dnorm) < 0)
324 return -1;
325 }
326
327 /* show after to catch 100% */
328 cursegm++;
329 if (totsegm < cursegm)
330 G_debug(1, "%d %d", totsegm, cursegm);
331
332 if (totsegm != 0) {
333 G_percent(cursegm, totsegm, 1);
334 }
335 /*
336 G_free_matrix(matrix);
337 G_free_ivector(indx);
338 G_free_vector(b);
339 */
340 G_free(data->points);
341 G_free(data);
342 }
343 return 1;
344}
345
346double smallest_segment(struct multtree *tree, int n_leafs)
347{
348 static int first_time = 1;
349 int ii;
350 static double minside;
351 double side;
352
353 if (tree == NULL)
354 return 0;
355 if (tree->data == NULL)
356 return 0;
357 if (tree->leafs != NULL) {
358 for (ii = 0; ii < n_leafs; ii++) {
359 side = smallest_segment(tree->leafs[ii], n_leafs);
360 if (first_time) {
361 minside = side;
362 first_time = 0;
363 }
364 if (side < minside)
365 minside = side;
366 }
367 }
368 else {
369 side = ((struct quaddata *)(tree->data))->xmax -
370 ((struct quaddata *)(tree->data))->x_orig;
371 return side;
372 }
373
374 return minside;
375}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#define NULL
Definition: ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition: dataquad.c:62
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double b
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition: ialloc.c:41
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:104
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition: qtree.c:194
double smallest_segment(struct multtree *tree, int n_leafs)
Definition: segmen2d.c:346
int IL_interp_segments_2d(struct interp_params *params, struct tree_info *info, struct multtree *tree, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot, int totsegm, off_t offset1, double dnorm)
Definition: segmen2d.c:47
check_points_fn * check_points
Definition: interpf.h:99
FILE * Tmp_fd_xx
Definition: interpf.h:93
FILE * Tmp_fd_xy
Definition: interpf.h:94
FILE * Tmp_fd_yy
Definition: interpf.h:94
grid_calc_fn * grid_calc
Definition: interpf.h:97
double x_orig
Definition: interpf.h:87
FILE * Tmp_fd_dx
Definition: interpf.h:92
double y_orig
Definition: interpf.h:87
FILE * Tmp_fd_z
Definition: interpf.h:92
FILE * Tmp_fd_dy
Definition: interpf.h:93
Definition: qtree.h:57
struct multtree ** leafs
Definition: qtree.h:59
struct quaddata * data
Definition: qtree.h:58
double ymax
Definition: dataquad.h:53
double y_orig
Definition: dataquad.h:51
double x_orig
Definition: dataquad.h:50
struct triple * points
Definition: dataquad.h:57
int n_points
Definition: dataquad.h:56
double xmax
Definition: dataquad.h:52
int n_cols
Definition: dataquad.h:55
int n_rows
Definition: dataquad.h:54
struct multtree * root
Definition: qtree.h:53
double z
Definition: dataquad.h:44
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43