GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
ressegm2d.c
Go to the documentation of this file.
1
2/*-
3 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
4 * University of Illinois
5 * US Army Construction Engineering Research Lab
6 * Copyright 1993, H. Mitasova (University of Illinois),
7 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
8 *
9 * modified by McCauley in August 1995
10 * modified by Mitasova in August 1995
11 *
12 * bug fixes by Jaro Hofierka in February 1999:
13 * line: 175,348 (*dnorm)
14 * 177,350 (points[m1].sm)
15 * 457,461 (})
16 *
17 * modified by Mitasova November 1999 (option for dnorm ind. tension)
18 *
19 */
20
21#include <stdio.h>
22#include <stdlib.h>
23#include <math.h>
24
25#include <grass/gis.h>
26#include <grass/raster.h>
27#include <grass/interpf.h>
28#include <grass/gmath.h>
29
30static int input_data(struct interp_params *,
31 int, int, struct fcell_triple *, int, int, int, int,
32 double, double, double);
33static int write_zeros(struct interp_params *, struct quaddata *, off_t);
34
35int IL_resample_interp_segments_2d(struct interp_params *params, struct BM *bitmask, /* bitmask */
36 double zmin, double zmax, /* min and max input z-values */
37 double *zminac, double *zmaxac, /* min and max interp. z-values */
38 double *gmin, double *gmax, /* min and max inperp. slope val. */
39 double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */
40 double *ertot, /* total interplating func. error */
41 off_t offset1, /* offset for temp file writing */
42 double *dnorm,
43 int overlap,
44 int inp_rows,
45 int inp_cols,
46 int fdsmooth,
47 int fdinp,
48 double ns_res,
49 double ew_res,
50 double inp_ns_res,
51 double inp_ew_res, int dtens)
52{
53
54 int i, j, k, l, m, m1, i1; /* loop coounters */
55 int cursegm = 0;
56 int new_comp = 0;
57 int n_rows, n_cols, inp_r, inp_c;
58 double x_or, y_or, xm, ym;
59 static int first = 1, new_first = 1;
60 double **matrix = NULL, **new_matrix = NULL, *b = NULL;
61 int *indx = NULL, *new_indx = NULL;
62 static struct fcell_triple *in_points = NULL; /* input points */
63 int inp_check_rows, inp_check_cols, /* total input rows/cols */
64 out_check_rows, out_check_cols; /* total output rows/cols */
65 int first_row, last_row; /* first and last input row of segment */
66 int first_col, last_col; /* first and last input col of segment */
67 int num, prev;
68 int div; /* number of divides */
69 int rem_out_row, rem_out_col; /* output rows/cols remainders */
70 int inp_seg_r, inp_seg_c, /* # of input rows/cols in segment */
71 out_seg_r, out_seg_c; /* # of output rows/cols in segment */
72 int ngstc, nszc /* first and last output col of the
73 * segment */
74 , ngstr, nszr; /* first and last output row of the
75 * segment */
76 int index; /* index for input data */
77 int c, r;
78 int overlap1;
79 int p_size;
80 struct quaddata *data;
81 double xmax, xmin, ymax, ymin;
82 int totsegm; /* total number of segments */
83 int total_points = 0;
84 struct triple triple; /* contains garbage */
85
86
87 xmin = params->x_orig;
88 ymin = params->y_orig;
89 xmax = xmin + ew_res * params->nsizc;
90 ymax = ymin + ns_res * params->nsizr;
91 prev = inp_rows * inp_cols;
92 if (prev <= params->kmax)
93 div = 1; /* no segmentation */
94
95 else { /* find the number of divides */
96 for (i = 2;; i++) {
97 c = inp_cols / i;
98 r = inp_rows / i;
99 num = c * r;
100 if (num < params->kmin) {
101 if (((params->kmin - num) > (prev + 1 - params->kmax)) &&
102 (prev + 1 < params->KMAX2)) {
103 div = i - 1;
104 break;
105 }
106 else {
107 div = i;
108 break;
109 }
110 }
111 if ((num > params->kmin) && (num + 1 < params->kmax)) {
112 div = i;
113 break;
114 }
115 prev = num;
116 }
117 }
118 out_seg_r = params->nsizr / div; /* output rows per segment */
119 out_seg_c = params->nsizc / div; /* output cols per segment */
120 inp_seg_r = inp_rows / div; /* input rows per segment */
121 inp_seg_c = inp_cols / div; /* input rows per segment */
122 rem_out_col = params->nsizc % div;
123 rem_out_row = params->nsizr % div;
124 overlap1 = min1(overlap, inp_seg_c - 1);
125 overlap1 = min1(overlap1, inp_seg_r - 1);
126 out_check_rows = 0;
127 out_check_cols = 0;
128 inp_check_rows = 0;
129 inp_check_cols = 0;
130
131 if (div == 1) {
132 p_size = inp_seg_c * inp_seg_r;
133 }
134 else {
135 p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
136 }
137 if (!in_points) {
138 if (!
139 (in_points =
140 (struct fcell_triple *)G_malloc(sizeof(struct fcell_triple) *
141 p_size * div))) {
142 fprintf(stderr, "Cannot allocate memory for in_points\n");
143 return -1;
144 }
145 }
146
147 *dnorm =
148 sqrt(((xmax - xmin) * (ymax -
149 ymin) * p_size) / (inp_rows * inp_cols));
150
151 if (dtens) {
152 params->fi = params->fi * (*dnorm) / 1000.;
153 fprintf(stderr, "dnorm = %f, rescaled tension = %f\n", *dnorm,
154 params->fi);
155 }
156
157 if (div == 1) { /* no segmentation */
158 totsegm = 1;
159 cursegm = 1;
160
161 input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
162 inp_cols, zmin, inp_ns_res, inp_ew_res);
163
164 x_or = 0.;
165 y_or = 0.;
166 xm = params->nsizc * ew_res;
167 ym = params->nsizr * ns_res;
168
169 data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym,
170 params->nsizr, params->nsizc,
171 0, params->KMAX2);
172 m1 = 0;
173 for (k = 1; k <= p_size; k++) {
174 if (!Rast_is_f_null_value(&(in_points[k - 1].z))) {
175 data->points[m1].x = in_points[k - 1].x / (*dnorm);
176 data->points[m1].y = in_points[k - 1].y / (*dnorm);
177 /* data->points[m1].z = (double) (in_points[k - 1].z) / (*dnorm); */
178 data->points[m1].z = (double)(in_points[k - 1].z);
179 data->points[m1].sm = in_points[k - 1].smooth;
180 m1++;
181 }
182 }
183 data->n_points = m1;
184 total_points = m1;
185 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
186 fprintf(stderr, "Cannot allocate memory for indx\n");
187 return -1;
188 }
189 if (!(matrix = G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
190 fprintf(stderr, "Cannot allocate memory for matrix\n");
191 return -1;
192 }
193 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
194 fprintf(stderr, "Cannot allocate memory for b\n");
195 return -1;
196 }
197
198 if (params->matrix_create(params, data->points, m1, matrix, indx) < 0)
199 return -1;
200 for (i = 0; i < m1; i++) {
201 b[i + 1] = data->points[i].z;
202 }
203 b[0] = 0.;
204 G_lubksb(matrix, m1 + 1, indx, b);
205
206 params->check_points(params, data, b, ertot, zmin, *dnorm, triple);
207
208 if (params->grid_calc(params, data, bitmask,
209 zmin, zmax, zminac, zmaxac, gmin, gmax,
210 c1min, c1max, c2min, c2max, ertot, b, offset1,
211 *dnorm) < 0) {
212 fprintf(stderr, "interpolation failed\n");
213 return -1;
214 }
215 else {
216 if (totsegm != 0) {
217 G_percent(cursegm, totsegm, 1);
218 }
219 /*
220 * if (b) G_free_vector(b); if (matrix) G_free_matrix(matrix); if
221 * (indx) G_free_ivector(indx);
222 */
223 fprintf(stderr, "dnorm in ressegm after grid before out= %f \n",
224 *dnorm);
225 return total_points;
226 }
227 }
228
229 out_seg_r = params->nsizr / div; /* output rows per segment */
230 out_seg_c = params->nsizc / div; /* output cols per segment */
231 inp_seg_r = inp_rows / div; /* input rows per segment */
232 inp_seg_c = inp_cols / div; /* input rows per segment */
233 rem_out_col = params->nsizc % div;
234 rem_out_row = params->nsizr % div;
235 overlap1 = min1(overlap, inp_seg_c - 1);
236 overlap1 = min1(overlap1, inp_seg_r - 1);
237 out_check_rows = 0;
238 out_check_cols = 0;
239 inp_check_rows = 0;
240 inp_check_cols = 0;
241
242 totsegm = div * div;
243
244 /* set up a segment */
245 for (i = 1; i <= div; i++) { /* input and output rows */
246 if (i <= div - rem_out_row)
247 n_rows = out_seg_r;
248 else
249 n_rows = out_seg_r + 1;
250 inp_r = inp_seg_r;
251 out_check_cols = 0;
252 inp_check_cols = 0;
253 ngstr = out_check_rows + 1; /* first output row of the segment */
254 nszr = ngstr + n_rows - 1; /* last output row of the segment */
255 y_or = (ngstr - 1) * ns_res; /* y origin of the segment */
256 /*
257 * Calculating input starting and ending rows and columns of this
258 * segment
259 */
260 first_row = (int)(y_or / inp_ns_res) + 1;
261 if (first_row > overlap1) {
262 first_row -= overlap1; /* middle */
263 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
264 if (last_row > inp_rows) {
265 first_row -= (last_row - inp_rows); /* bottom */
266 last_row = inp_rows;
267 }
268 }
269 else {
270 first_row = 1; /* top */
271 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
272 }
273 if ((last_row > inp_rows) || (first_row < 1)) {
274 fprintf(stderr, "Row overlap too large!\n");
275 return -1;
276 }
277 input_data(params, first_row, last_row, in_points, fdsmooth, fdinp,
278 inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);
279
280 for (j = 1; j <= div; j++) { /* input and output cols */
281 if (j <= div - rem_out_col)
282 n_cols = out_seg_c;
283 else
284 n_cols = out_seg_c + 1;
285 inp_c = inp_seg_c;
286
287 ngstc = out_check_cols + 1; /* first output col of the segment */
288 nszc = ngstc + n_cols - 1; /* last output col of the segment */
289 x_or = (ngstc - 1) * ew_res; /* x origin of the segment */
290
291 first_col = (int)(x_or / inp_ew_res) + 1;
292 if (first_col > overlap1) {
293 first_col -= overlap1; /* middle */
294 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
295 if (last_col > inp_cols) {
296 first_col -= (last_col - inp_cols); /* right */
297 last_col = inp_cols;
298 }
299 }
300 else {
301 first_col = 1; /* left */
302 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
303 }
304 if ((last_col > inp_cols) || (first_col < 1)) {
305 fprintf(stderr, "Column overlap too large!\n");
306 return -1;
307 }
308 m = 0;
309 /* Getting points for interpolation (translated) */
310
311 xm = nszc * ew_res;
312 ym = nszr * ns_res;
313 data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym,
314 nszr - ngstr + 1,
315 nszc - ngstc + 1, 0,
316 params->KMAX2);
317 new_comp = 0;
318
319 for (k = 0; k <= last_row - first_row; k++) {
320 for (l = first_col - 1; l < last_col; l++) {
321 index = k * inp_cols + l;
322 if (!Rast_is_f_null_value(&(in_points[index].z))) {
323 /* if the point is inside the segment (not overlapping) */
324 if ((in_points[index].x - x_or >= 0) &&
325 (in_points[index].y - y_or >= 0) &&
326 ((nszc - 1) * ew_res - in_points[index].x >= 0) &&
327 ((nszr - 1) * ns_res - in_points[index].y >= 0))
328 total_points += 1;
329 data->points[m].x =
330 (in_points[index].x - x_or) / (*dnorm);
331 data->points[m].y =
332 (in_points[index].y - y_or) / (*dnorm);
333 /* data->points[m].z = (double) (in_points[index].z) / (*dnorm); */
334 data->points[m].z = (double)(in_points[index].z);
335 data->points[m].sm = in_points[index].smooth;
336 m++;
337 }
338 else
339 new_comp = 1;
340
341 /* fprintf(stderr,"%f,%f,%f
342 zmin=%f\n",in_points[index].x,in_points[index].y,in_points[index].z,zmin);
343 */
344 }
345 }
346 /* fprintf (stdout,"m,index:%di,%d\n",m,index); */
347 if (m <= params->KMAX2)
348 data->n_points = m;
349 else
350 data->n_points = params->KMAX2;
351 out_check_cols += n_cols;
352 inp_check_cols += inp_c;
353 cursegm = (i - 1) * div + j - 1;
354
355 /* show before to catch 0% */
356 if (totsegm != 0) {
357 G_percent(cursegm, totsegm, 1);
358 }
359 if (m == 0) {
360 /*
361 * fprintf(stderr,"Warning: segment with zero points encountered,
362 * insrease overlap\n");
363 */
364 write_zeros(params, data, offset1);
365 }
366 else {
367 if (new_comp) {
368 if (new_first) {
369 new_first = 0;
370 if (!b) {
371 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
372 fprintf(stderr,
373 "Cannot allocate memory for b\n");
374 return -1;
375 }
376 }
377 if (!(new_indx = G_alloc_ivector(params->KMAX2 + 1))) {
378 fprintf(stderr,
379 "Cannot allocate memory for new_indx\n");
380 return -1;
381 }
382 if (!
383 (new_matrix =
384 G_alloc_matrix(params->KMAX2 + 1,
385 params->KMAX2 + 1))) {
386 fprintf(stderr,
387 "Cannot allocate memory for new_matrix\n");
388 return -1;
389 }
390 } /*new_first */
391 if (params->
392 matrix_create(params, data->points, data->n_points,
393 new_matrix, new_indx) < 0)
394 return -1;
395
396 for (i1 = 0; i1 < m; i1++) {
397 b[i1 + 1] = data->points[i1].z;
398 }
399 b[0] = 0.;
400 G_lubksb(new_matrix, data->n_points + 1, new_indx, b);
401
402 params->check_points(params, data, b, ertot, zmin,
403 *dnorm, triple);
404
405 if (params->grid_calc(params, data, bitmask,
406 zmin, zmax, zminac, zmaxac, gmin,
407 gmax, c1min, c1max, c2min, c2max,
408 ertot, b, offset1, *dnorm) < 0) {
409
410 fprintf(stderr, "interpolate() failed\n");
411 return -1;
412 }
413 } /*new_comp */
414 else {
415 if (first) {
416 first = 0;
417 if (!b) {
418 if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
419 fprintf(stderr,
420 "Cannot allocate memory for b\n");
421 return -1;
422 }
423 }
424 if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
425 fprintf(stderr,
426 "Cannot allocate memory for indx\n");
427 return -1;
428 }
429 if (!
430 (matrix =
431 G_alloc_matrix(params->KMAX2 + 1,
432 params->KMAX2 + 1))) {
433 fprintf(stderr,
434 "Cannot allocate memory for matrix\n");
435 return -1;
436 }
437 } /* first */
438 if (params->
439 matrix_create(params, data->points, data->n_points,
440 matrix, indx) < 0)
441 return -1;
442 /* } here it was bug */
443 for (i1 = 0; i1 < m; i1++)
444 b[i1 + 1] = data->points[i1].z;
445 b[0] = 0.;
446 G_lubksb(matrix, data->n_points + 1, indx, b);
447
448 params->check_points(params, data, b, ertot, zmin,
449 *dnorm, triple);
450
451 if (params->grid_calc(params, data, bitmask,
452 zmin, zmax, zminac, zmaxac, gmin,
453 gmax, c1min, c1max, c2min, c2max,
454 ertot, b, offset1, *dnorm) < 0) {
455
456 fprintf(stderr, "interpolate() failed\n");
457 return -1;
458 }
459 }
460 }
461 if (data) {
462 G_free(data->points);
463 G_free(data);
464 }
465 /*
466 * cursegm++;
467 */
468 }
469
470 inp_check_rows += inp_r;
471 out_check_rows += n_rows;
472 }
473
474 /* run one last time after the loop is done to catch 100% */
475 if (totsegm != 0)
476 G_percent(1, 1, 1); /* cursegm doesn't get to totsegm so we force 100% */
477
478 /*
479 * if (b) G_free_vector(b); if (indx) G_free_ivector(indx); if (matrix)
480 * G_free_matrix(matrix);
481 */
482 fprintf(stderr, "dnorm in ressegm after grid before out2= %f \n", *dnorm);
483 return total_points;
484}
485
486/* input of data for interpolation and smoothing parameters */
487
488static int input_data(struct interp_params *params,
489 int first_row, int last_row,
490 struct fcell_triple *points,
491 int fdsmooth, int fdinp,
492 int inp_rows, int inp_cols,
493 double zmin, double inp_ns_res, double inp_ew_res)
494{
495 double x, y, sm; /* input data and smoothing */
496 int m1, m2; /* loop counters */
497 static FCELL *cellinp = NULL; /* cell buffer for input data */
498 static FCELL *cellsmooth = NULL; /* cell buffer for smoothing */
499
500
501 if (!cellinp)
502 cellinp = Rast_allocate_f_buf();
503 if (!cellsmooth)
504 cellsmooth = Rast_allocate_f_buf();
505
506 for (m1 = 0; m1 <= last_row - first_row; m1++) {
507 Rast_get_f_row(fdinp, cellinp, inp_rows - m1 - first_row);
508 if (fdsmooth >= 0)
509 Rast_get_f_row(fdsmooth, cellsmooth, inp_rows - m1 - first_row);
510
511 y = params->y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res;
512 for (m2 = 0; m2 < inp_cols; m2++) {
513 x = params->x_orig + (m2 + 0.5) * inp_ew_res;
514 /*
515 * z = cellinp[m2]*params->zmult;
516 */
517 if (fdsmooth >= 0)
518 sm = (double)cellsmooth[m2];
519 else
520 sm = 0.01;
521
522 points[m1 * inp_cols + m2].x = x - params->x_orig;
523 points[m1 * inp_cols + m2].y = y - params->y_orig;
524 if (!Rast_is_f_null_value(cellinp + m2)) {
525 points[m1 * inp_cols + m2].z =
526 cellinp[m2] * params->zmult - zmin;
527 }
528 else {
529 Rast_set_f_null_value(&(points[m1 * inp_cols + m2].z), 1);
530 }
531
532 /* fprintf (stdout,"sm: %f\n",sm); */
533
534 points[m1 * inp_cols + m2].smooth = sm;
535 }
536 }
537 return 1;
538}
539
540static int write_zeros(struct interp_params *params,
541 struct quaddata *data, /* given segment */
542 off_t offset1 /* offset for temp file writing */
543 )
544{
545
546 /*
547 * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
548 * c
549 */
550 double x_or = data->x_orig;
551 double y_or = data->y_orig;
552 int n_rows = data->n_rows;
553 int n_cols = data->n_cols;
554 int cond1, cond2;
555 int k, l;
556 int ngstc, nszc, ngstr, nszr;
557 off_t offset, offset2;
558 double ns_res, ew_res;
559
560 ns_res = (((struct quaddata *)(data))->ymax -
561 ((struct quaddata *)(data))->y_orig) / data->n_rows;
562 ew_res = (((struct quaddata *)(data))->xmax -
563 ((struct quaddata *)(data))->x_orig) / data->n_cols;
564
565 cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
566 (params->adxy != NULL));
567 cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);
568
569 ngstc = (int)(x_or / ew_res + 0.5) + 1;
570 nszc = ngstc + n_cols - 1;
571 ngstr = (int)(y_or / ns_res + 0.5) + 1;
572 nszr = ngstr + n_rows - 1;
573
574 for (k = ngstr; k <= nszr; k++) {
575 offset = offset1 * (k - 1); /* rows offset */
576 for (l = ngstc; l <= nszc; l++) {
577 /*
578 * params->az[l] = 0.;
579 */
580 Rast_set_d_null_value(params->az + l, 1);
581 if (cond1) {
582 /*
583 * params->adx[l] = (FCELL)0.; params->ady[l] = (FCELL)0.;
584 */
585 Rast_set_d_null_value(params->adx + l, 1);
586 Rast_set_d_null_value(params->ady + l, 1);
587 if (cond2) {
588 Rast_set_d_null_value(params->adxx + l, 1);
589 Rast_set_d_null_value(params->adyy + l, 1);
590 Rast_set_d_null_value(params->adxy + l, 1);
591 /*
592 * params->adxx[l] = (FCELL)0.; params->adyy[l] = (FCELL)0.;
593 * params->adxy[l] = (FCELL)0.;
594 */
595 }
596 }
597 }
598 offset2 = (offset + ngstc - 1) * sizeof(FCELL);
599 if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
600 return -1;
601 }
602 return 1;
603}
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
double b
double l
double r
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition: ialloc.c:41
int min1(int, int)
Definition: minmax.c:18
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 IL_resample_interp_segments_2d(struct interp_params *params, 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, off_t offset1, double *dnorm, int overlap, int inp_rows, int inp_cols, int fdsmooth, int fdinp, double ns_res, double ew_res, double inp_ns_res, double inp_ew_res, int dtens)
Definition: ressegm2d.c:35
double x
Definition: interpf.h:18
double smooth
Definition: interpf.h:21
double y
Definition: interpf.h:19
check_points_fn * check_points
Definition: interpf.h:99
double zmult
Definition: interpf.h:70
DCELL * az
Definition: interpf.h:78
grid_calc_fn * grid_calc
Definition: interpf.h:97
double fi
Definition: interpf.h:80
double x_orig
Definition: interpf.h:87
DCELL * adxy
Definition: interpf.h:79
DCELL * adyy
Definition: interpf.h:79
DCELL * adx
Definition: interpf.h:78
double y_orig
Definition: interpf.h:87
DCELL * ady
Definition: interpf.h:78
wr_temp_fn * wr_temp
Definition: interpf.h:103
DCELL * adxx
Definition: interpf.h:79
matrix_create_fn * matrix_create
Definition: interpf.h:98
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
double z
Definition: dataquad.h:44
double sm
Definition: dataquad.h:45
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43
#define x