GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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#include <stdio.h>
18#include <stdlib.h>
19#include <math.h>
20#include <grass/gis.h>
21#include <grass/glocale.h>
22#include <grass/interpf.h>
23#include <grass/gmath.h>
24
25/*!
26 * Interpolate recursively a tree of segments
27 *
28 * Recursively processes each segment in a tree by:
29 * - finding points from neighbouring segments so that the total number of
30 * points is between KMIN and KMAX2 by calling tree function MT_get_region().
31 * - creating and solving the system of linear equations using these points
32 * and interp() by calling matrix_create() and G_ludcmp().
33 * - checking the interpolating function values at points by calling
34 * check_points().
35 * - computing grid for this segment using points and interp() by calling
36 * grid_calc().
37 *
38 * \todo
39 * Isn't this in fact the updated version of the function
40 * (IL_interp_segments_new_2d)? The function IL_interp_segments_new_2d has the
41 * following, better behavior: The difference between this function and
42 * IL_interp_segments_2d() is making sure that additional points are taken from
43 * all directions, i.e. it finds equal number of points from neighboring
44 * segments in each of 8 neighborhoods.
45 */
47 struct interp_params *params,
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, 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 = (((struct quaddata *)(info->root->data))->xmax -
87 ((struct quaddata *)(info->root->data))->x_orig) /
88 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], bitmask, zmin,
97 zmax, zminac, zmaxac, gmin, gmax, c1min,
98 c1max, c2min, c2max, ertot, totsegm, offset1,
99 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
115 this case 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 = params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
122 /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf,
123 * DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
124
125 data = (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
126 xmx + distx, ymx + disty, 0, 0,
127 0, params->KMAX2);
128 npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
129
130 while ((npt < MINPTS) || (npt > params->KMAX2)) {
131 if (i >= 70) {
132 G_warning(
133 _("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 temp1 = distxp;
144 distxp = distx;
145 distx = distxp - fabs(distx - temp1) * 0.5;
146 temp2 = distyp;
147 distyp = disty;
148 disty = distyp - fabs(disty - temp2) * 0.5;
149 /* decrease by 50% of a previous change in window */
150 }
151 else {
152 temp1 = distyp;
153 distyp = disty;
154 temp2 = distxp;
155 distxp = distx;
156 if (MAXENC) {
157 disty = fabs(disty - temp1) * 0.5 + distyp;
158 distx = fabs(distx - temp2) * 0.5 + distxp;
159 }
160 else {
161 distx += distx;
162 disty += disty;
163 }
164 /* decrease by 50% of extra distance */
165 }
166 data->x_orig = xmn - distx; /* update window */
167 data->y_orig = ymn - disty;
168 data->xmax = xmx + distx;
169 data->ymax = ymx + disty;
170 data->n_points = 0;
171 npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
172 }
173
174 if (totsegm != 0) {
175 G_percent(cursegm, totsegm, 1);
176 }
177 data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
178 data->n_cols = ((struct quaddata *)(tree->data))->n_cols;
179
180 /* for printing out overlapping segments */
181 ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
182 ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
183 ((struct quaddata *)(tree->data))->xmax = xmx + distx;
184 ((struct quaddata *)(tree->data))->ymax = ymx + disty;
185
186 data->x_orig = xmn;
187 data->y_orig = ymn;
188 data->xmax = xmx;
189 data->ymax = ymx;
190
191 if (!matrix) {
192 if (!(matrix =
193 G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
194 G_warning(_("Out of memory"));
195 return -1;
196 }
197 }
198 if (!indx) {
199 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
200 G_warning(_("Out of memory"));
201 return -1;
202 }
203 }
204 if (!b) {
205 if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
206 G_warning(_("Out of memory"));
207 return -1;
208 }
209 }
210 /* allocate memory for CV points only if cv is performed */
211 if (params->cv) {
212 if (!(point = (struct triple *)G_malloc(sizeof(struct triple) *
213 data->n_points))) {
214 G_warning(_("Out of memory"));
215 return -1;
216 }
217 }
218
219 /*normalize the data so that the side of average segment is about 1m */
220 /* put data_points into point only if CV is performed */
221
222 for (i = 0; i < data->n_points; i++) {
223 data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
224 data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
225 if (params->cv) {
226 point[i].x = data->points[i].x; /*cv stuff */
227 point[i].y = data->points[i].y; /*cv stuff */
228 point[i].z = data->points[i].z; /*cv stuff */
229 }
230
231 /* commented out by Helena january 1997 as this is not necessary
232 although it may be useful to put normalization of z back?
233 data->points[i].z = data->points[i].z / dnorm;
234 this made smoothing self-adjusting based on dnorm
235 if (params->rsm < 0.) data->points[i].sm = data->points[i].sm /
236 dnorm;
237 */
238 }
239
240 /* cv stuff */
241 if (params->cv)
242 m_skip = data->n_points;
243 else
244 m_skip = 1;
245
246 /* remove after cleanup - this is just for testing */
247 skip_point.x = 0.;
248 skip_point.y = 0.;
249 skip_point.z = 0.;
250
251 /*** TODO: parallelize this loop instead of the LU solver! ***/
252 for (skip_index = 0; skip_index < m_skip; skip_index++) {
253 if (params->cv) {
254 segtest = 0;
255 j = 0;
256 xx =
257 point[skip_index].x * dnorm + data->x_orig + params->x_orig;
258 yy =
259 point[skip_index].y * dnorm + data->y_orig + params->y_orig;
260 /* zz = point[skip_index].z; */
261 if (xx >= data->x_orig + params->x_orig &&
262 xx <= data->xmax + params->x_orig &&
263 yy >= data->y_orig + params->y_orig &&
264 yy <= data->ymax + params->y_orig) {
265 segtest = 1;
266 skip_point.x = point[skip_index].x;
267 skip_point.y = point[skip_index].y;
268 skip_point.z = point[skip_index].z;
269 for (k = 0; k < m_skip; k++) {
270 if (k != skip_index && params->cv) {
271 data->points[j].x = point[k].x;
272 data->points[j].y = point[k].y;
273 data->points[j].z = point[k].z;
274 j++;
275 }
276 }
277 } /* segment area test */
278 }
279 if (!params->cv) {
280 if (params->matrix_create(params, data->points, data->n_points,
281 matrix, indx) < 0)
282 return -1;
283 }
284 else if (segtest == 1) {
285 if (params->matrix_create(params, data->points,
286 data->n_points - 1, matrix, indx) < 0)
287 return -1;
288 }
289 if (!params->cv) {
290 for (i = 0; i < data->n_points; i++)
291 b[i + 1] = data->points[i].z;
292 b[0] = 0.;
293 G_lubksb(matrix, data->n_points + 1, indx, b);
294 /* put here condition to skip error if not needed */
295 params->check_points(params, data, b, ertot, zmin, dnorm,
296 skip_point);
297 }
298 else if (segtest == 1) {
299 for (i = 0; i < data->n_points - 1; i++)
300 b[i + 1] = data->points[i].z;
301 b[0] = 0.;
302 G_lubksb(matrix, data->n_points, indx, b);
303 params->check_points(params, data, b, ertot, zmin, dnorm,
304 skip_point);
305 }
306 } /*end of cv loop */
307
308 if (!params->cv)
309 if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
310 (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
311 (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {
312
313 if (params->grid_calc(params, data, bitmask, zmin, zmax, zminac,
314 zmaxac, gmin, gmax, c1min, c1max, c2min,
315 c2max, ertot, b, offset1, dnorm) < 0)
316 return -1;
317 }
318
319 /* show after to catch 100% */
320 cursegm++;
321 if (totsegm < cursegm)
322 G_debug(1, "%d %d", totsegm, cursegm);
323
324 if (totsegm != 0) {
325 G_percent(cursegm, totsegm, 1);
326 }
327 /*
328 G_free_matrix(matrix);
329 G_free_ivector(indx);
330 G_free_vector(b);
331 */
332 G_free(data->points);
333 G_free(data);
334 }
335 return 1;
336}
337
338double smallest_segment(struct multtree *tree, int n_leafs)
339{
340 static int first_time = 1;
341 int ii;
342 static double minside;
343 double side;
344
345 if (tree == NULL)
346 return 0;
347 if (tree->data == NULL)
348 return 0;
349 if (tree->leafs != NULL) {
350 for (ii = 0; ii < n_leafs; ii++) {
351 side = smallest_segment(tree->leafs[ii], n_leafs);
352 if (first_time) {
353 minside = side;
354 first_time = 0;
355 }
356 if (side < minside)
357 minside = side;
358 }
359 }
360 else {
361 side = ((struct quaddata *)(tree->data))->xmax -
362 ((struct quaddata *)(tree->data))->x_orig;
363 return side;
364 }
365
366 return minside;
367}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define NULL
Definition ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:57
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:39
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:59
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double b
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition ialloc.c:38
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition lu.c:103
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition qtree.c:186
double smallest_segment(struct multtree *tree, int n_leafs)
Definition segmen2d.c:338
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:46
check_points_fn * check_points
Definition interpf.h:120
FILE * Tmp_fd_xx
Definition interpf.h:111
FILE * Tmp_fd_xy
Definition interpf.h:111
FILE * Tmp_fd_yy
Definition interpf.h:111
grid_calc_fn * grid_calc
Definition interpf.h:116
double x_orig
Definition interpf.h:101
FILE * Tmp_fd_dx
Definition interpf.h:111
double y_orig
Definition interpf.h:101
FILE * Tmp_fd_z
Definition interpf.h:111
FILE * Tmp_fd_dy
Definition interpf.h:111
matrix_create_fn * matrix_create
Definition interpf.h:118
struct multtree ** leafs
Definition qtree.h:54
struct quaddata * data
Definition qtree.h:53
double ymax
Definition dataquad.h:49
double y_orig
Definition dataquad.h:47
double x_orig
Definition dataquad.h:46
struct triple * points
Definition dataquad.h:53
int n_points
Definition dataquad.h:52
double xmax
Definition dataquad.h:48
int n_cols
Definition dataquad.h:51
int n_rows
Definition dataquad.h:50
struct multtree * root
Definition qtree.h:49
double z
Definition dataquad.h:41
double x
Definition dataquad.h:39
double y
Definition dataquad.h:40