GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
gvl_calc.c
Go to the documentation of this file.
1/*!
2 \file lib/ogsf/gvl_calc.c
3
4 \brief OGSF library - loading and manipulating volumes (lower level
5 functions)
6
7 GRASS OpenGL gsurf OGSF Library
8
9 (C) 1999-2008 by the GRASS Development Team
10
11 This program is free software under the
12 GNU General Public License (>=v2).
13 Read the file COPYING that comes with GRASS
14 for details.
15
16 \author Tomas Paudits (February 2004)
17 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
18 */
19
20#include <math.h>
21
22#include <grass/gis.h>
23#include <grass/ogsf.h>
24
25#include "rgbpack.h"
26#include "mc33_table.h"
27
28/*!
29 \brief memory buffer for writing
30 */
31#define BUFFER_SIZE 1000000
32
33/* USEFUL MACROS */
34
35/* interp. */
36#define LINTERP(d, a, b) (a + d * (b - a))
37#define TINTERP(d, v) \
38 ((v[0] * (1. - d[0]) * (1. - d[1]) * (1. - d[2])) + \
39 (v[1] * d[0] * (1. - d[1]) * (1. - d[2])) + \
40 (v[2] * d[0] * d[1] * (1. - d[2])) + \
41 (v[3] * (1. - d[0]) * d[1] * (1. - d[2])) + \
42 (v[4] * (1. - d[0]) * (1. - d[1]) * d[2]) + \
43 (v[5] * d[0] * (1. - d[1]) * d[2]) + (v[6] * d[0] * d[1] * d[2]) + \
44 (v[7] * (1. - d[0]) * d[1] * d[2]))
45
46#define FOR_VAR i_for
47#define FOR_0_TO_N(n, cmd) \
48 { \
49 int FOR_VAR; \
50 for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) { \
51 cmd; \
52 } \
53 }
54
55/*!
56 \brief writing and reading isosurface data
57 */
58#define WRITE(c) gvl_write_char(dbuff->ndx_new++, &(dbuff->new), c)
59#define READ() gvl_read_char(dbuff->ndx_old++, dbuff->old)
60#define SKIP(n) dbuff->ndx_old = dbuff->ndx_old + n
61
62/*!
63 \brief check and set data descriptor
64 */
65#define IS_IN_DATA(att) ((isosurf->data_desc >> att) & 1)
66#define SET_IN_DATA(att) isosurf->data_desc = (isosurf->data_desc | (1 << att))
67
68typedef struct {
69 unsigned char *old;
70 unsigned char *new;
71 int ndx_old;
72 int ndx_new;
73 int num_zero;
74} data_buffer;
75
76int mc33_process_cube(int c_ndx, float *v);
77
78/* global variables */
80double ResX, ResY, ResZ;
81
82/************************************************************************/
83/* ISOSURFACES */
84
85/*!
86 \brief Write cube index
87
88 \param ndx
89 \param dbuff
90 */
91void iso_w_cndx(int ndx, data_buffer *dbuff)
92{
93 /* cube don't contains polys */
94 if (ndx == -1) {
95 if (dbuff->num_zero == 0) {
96 WRITE(0);
97 dbuff->num_zero++;
98 }
99 else if (dbuff->num_zero == 254) {
100 WRITE(dbuff->num_zero + 1);
101 dbuff->num_zero = 0;
102 }
103 else {
104 dbuff->num_zero++;
105 }
106 }
107 else { /* isosurface cube */
108 if (dbuff->num_zero == 0) {
109 WRITE((ndx / 256) + 1);
110 WRITE(ndx % 256);
111 }
112 else {
113 WRITE(dbuff->num_zero);
114 dbuff->num_zero = 0;
115 WRITE((ndx / 256) + 1);
116 WRITE(ndx % 256);
117 }
118 }
119}
120
121/*!
122 \brief Read cube index
123
124 \param dbuff
125 */
126int iso_r_cndx(data_buffer *dbuff)
127{
128 int ndx, ndx2;
129
130 if (dbuff->num_zero != 0) {
131 dbuff->num_zero--;
132 ndx = -1;
133 }
134 else {
135 WRITE(ndx = READ());
136 if (ndx == 0) {
137 WRITE(dbuff->num_zero = READ());
138 dbuff->num_zero--;
139 ndx = -1;
140 }
141 else {
142 WRITE(ndx2 = READ());
143 ndx = (ndx - 1) * 256 + ndx2;
144 }
145 }
146
147 return ndx;
148}
149
150/*!
151 \brief Get value from data input
152
153 \param isosurf
154 \param desc
155 \param x,y,z
156 \param[out] value
157
158 \return 0
159 \return ?
160 */
161int iso_get_cube_value(geovol_isosurf *isosurf, int desc, int x, int y, int z,
162 float *v)
163{
164 double d;
165 geovol_file *vf;
166 int type, ret = 1;
167
168 /* get volume file from attribute handle */
169 vf = gvl_file_get_volfile(isosurf->att[desc].hfile);
170 type = gvl_file_get_data_type(vf);
171
172 /* get value from volume file */
173 if (type == VOL_DTYPE_FLOAT) {
174 gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
175 (int)(z * ResZ), v);
176 }
177 else if (type == VOL_DTYPE_DOUBLE) {
178 gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
179 (int)(z * ResZ), &d);
180 *v = (float)d;
181 }
182 else {
183 return 0;
184 }
185
186 /* null check */
187 if (gvl_file_is_null_value(vf, v))
188 ret = 0;
189
190 /* adjust data */
191 switch (desc) {
192 case (ATT_TOPO):
193 *v = (*v) - isosurf->att[desc].constant;
194 break;
195 case (ATT_MASK):
196 if (isosurf->att[desc].constant)
197 ret = !ret;
198 break;
199 }
200
201 return ret;
202}
203
204/*!
205 \brief Get volume file values range
206
207 \param isosurf
208 \param desc
209 \param[out] min
210 \param[out] max
211 */
212void iso_get_range(geovol_isosurf *isosurf, int desc, double *min, double *max)
213{
214 gvl_file_get_min_max(gvl_file_get_volfile(isosurf->att[desc].hfile), min,
215 max);
216}
217
218/*!
219 \brief Read values for cube
220
221 \param isosurf
222 \param desc
223 \param x,y,z
224 \param[out] v
225
226 \return
227 */
228int iso_get_cube_values(geovol_isosurf *isosurf, int desc, int x, int y, int z,
229 float *v)
230{
231 int p, ret = 1;
232
233 for (p = 0; p < 8; ++p) {
234 if (iso_get_cube_value(isosurf, desc, x + ((p ^ (p >> 1)) & 1),
235 y + ((p >> 1) & 1), z + ((p >> 2) & 1),
236 &v[p]) == 0) {
237 ret = 0;
238 }
239 }
240
241 return ret;
242}
243
244/*!
245 \brief Calculate cube grads
246
247 \param isosurf
248 \param x,y,z
249 \param grad
250 */
251void iso_get_cube_grads(geovol_isosurf *isosurf, int x, int y, int z,
252 float (*grad)[3])
253{
254 float v[3];
255 int i, j, k, p;
256
257 for (p = 0; p < 8; ++p) {
258 i = x + ((p ^ (p >> 1)) & 1);
259 j = y + ((p >> 1) & 1);
260 k = z + ((p >> 2) & 1);
261
262 /* x */
263 if (i == 0) {
264 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
265 iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
266 grad[p][0] = v[2] - v[1];
267 }
268 else {
269 if (i == (Cols - 1)) {
270 iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
271 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
272 grad[p][0] = v[1] - v[0];
273 }
274 else {
275 iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
276 iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
277 grad[p][0] = (v[2] - v[0]) / 2;
278 }
279 }
280
281 /* y */
282 if (j == 0) {
283 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
284 iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
285 grad[p][1] = v[2] - v[1];
286 }
287 else {
288 if (j == (Rows - 1)) {
289 iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
290 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
291 grad[p][1] = v[1] - v[0];
292 }
293 else {
294 iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
295 iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
296 grad[p][1] = (v[2] - v[0]) / 2;
297 }
298 }
299
300 /* z */
301 if (k == 0) {
302 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
303 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
304 grad[p][2] = v[2] - v[1];
305 }
306 else {
307 if (k == (Depths - 1)) {
308 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
309 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
310 grad[p][2] = v[1] - v[0];
311 }
312 else {
313 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
314 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
315 grad[p][2] = (v[2] - v[0]) / 2;
316 }
317 }
318 }
319}
320
321/*!
322 \brief Process cube
323
324 \param isosurf
325 \param x,y,z
326 \param dbuff
327 */
328void iso_calc_cube(geovol_isosurf *isosurf, int x, int y, int z,
329 data_buffer *dbuff)
330{
331 int i, c_ndx;
332 int crnt, v1, v2, c;
333 float val[MAX_ATTS][8], grad[8][3];
334 float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
335 double min, max;
336
337 if (isosurf->att[ATT_TOPO].changed) {
338 /* read topo values, if there are NULL values then return */
339 if (!iso_get_cube_values(isosurf, ATT_TOPO, x, y, z, val[ATT_TOPO])) {
340 iso_w_cndx(-1, dbuff);
341 return;
342 }
343
344 /* mask */
345 if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
346 if (!iso_get_cube_values(isosurf, ATT_MASK, x, y, z,
347 val[ATT_MASK])) {
348 iso_w_cndx(-1, dbuff);
349 return;
350 }
351 }
352
353 /* index to precalculated table */
354 c_ndx = 0;
355 for (i = 0; i < 8; i++) {
356 if (val[ATT_TOPO][i] > 0)
357 c_ndx |= 1 << i;
358 }
359 c_ndx = mc33_process_cube(c_ndx, val[ATT_TOPO]);
360
361 iso_w_cndx(c_ndx, dbuff);
362
363 if (c_ndx == -1)
364 return;
365
366 /* calc cube grads */
367 iso_get_cube_grads(isosurf, x, y, z, grad);
368 }
369 else {
370 /* read cube index */
371 if ((c_ndx = iso_r_cndx(dbuff)) == -1)
372 return;
373 }
374
375 /* get color values */
376 if (isosurf->att[ATT_COLOR].changed &&
377 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
378 iso_get_cube_values(isosurf, ATT_COLOR, x, y, z, val[ATT_COLOR]);
379 }
380
381 /* get transparency values */
382 if (isosurf->att[ATT_TRANSP].changed &&
383 isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
384 iso_get_cube_values(isosurf, ATT_TRANSP, x, y, z, val[ATT_TRANSP]);
385 }
386
387 /* get shine values */
388 if (isosurf->att[ATT_SHINE].changed &&
389 isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
390 iso_get_cube_values(isosurf, ATT_SHINE, x, y, z, val[ATT_SHINE]);
391 }
392
393 /* get emit values */
394 if (isosurf->att[ATT_EMIT].changed &&
395 isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
396 iso_get_cube_values(isosurf, ATT_EMIT, x, y, z, val[ATT_EMIT]);
397 }
398
399 FOR_0_TO_N(3, d_sum[FOR_VAR] = 0.; n_sum[FOR_VAR] = 0.);
400
401 /* loop in edges */
402 for (i = 0; i < cell_table[c_ndx].nedges; i++) {
403 /* get edge number */
404 crnt = cell_table[c_ndx].edges[i];
405
406 /* set topo */
407 if (isosurf->att[ATT_TOPO].changed) {
408 /* interior vertex */
409 if (crnt == 12) {
410 FOR_0_TO_N(3, WRITE((d3[FOR_VAR] =
411 d_sum[FOR_VAR] /
412 ((float)(cell_table[c_ndx].nedges))) *
413 255));
414 GS_v3norm(n_sum);
415 FOR_0_TO_N(3, WRITE((n_sum[FOR_VAR] /
416 ((float)(cell_table[c_ndx].nedges)) +
417 1.) *
418 127));
419 /* edge vertex */
420 }
421 else {
422 /* set edges verts */
423 v1 = edge_vert[crnt][0];
424 v2 = edge_vert[crnt][1];
425
426 /* calc intersection point - edge and isosurf */
427 d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] - val[ATT_TOPO][v2]);
428
429 d_sum[edge_vert_pos[crnt][0]] += d;
430 d_sum[edge_vert_pos[crnt][1]] += edge_vert_pos[crnt][2];
431 d_sum[edge_vert_pos[crnt][3]] += edge_vert_pos[crnt][4];
432
433 WRITE(d * 255);
434
435 /* set normal for intersect. point */
436 FOR_0_TO_N(3, n[FOR_VAR] = LINTERP(d, grad[v1][FOR_VAR],
437 grad[v2][FOR_VAR]));
438 GS_v3norm(n);
439 FOR_0_TO_N(3, n_sum[FOR_VAR] += n[FOR_VAR]);
440 FOR_0_TO_N(3, WRITE((n[FOR_VAR] + 1.) * 127));
441 }
442 }
443 else {
444 /* read x,y,z of intersection point in cube coords */
445 if (crnt == 12) {
446 WRITE(c = READ());
447 d3[0] = ((float)c) / 255.0;
448 WRITE(c = READ());
449 d3[1] = ((float)c) / 255.0;
450 WRITE(c = READ());
451 d3[2] = ((float)c) / 255.0;
452 }
453 else {
454 /* set edges verts */
455 v1 = edge_vert[crnt][0];
456 v2 = edge_vert[crnt][1];
457
458 WRITE(c = READ());
459 d = ((float)c) / 255.0;
460 }
461
462 /* set normals */
463 FOR_0_TO_N(3, WRITE(READ()));
464 }
465
466 /* set color */
467 if (isosurf->att[ATT_COLOR].changed &&
468 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
469 if (crnt == 12) {
470 tv = TINTERP(d3, val[ATT_COLOR]);
471 }
472 else {
473 tv = LINTERP(d, val[ATT_COLOR][v1], val[ATT_COLOR][v2]);
474 }
475
476 c = Gvl_get_color_for_value(isosurf->att[ATT_COLOR].att_data, &tv);
477
478 WRITE(c & RED_MASK);
479 WRITE((c & GRN_MASK) >> 8);
480 WRITE((c & BLU_MASK) >> 16);
481
482 if (IS_IN_DATA(ATT_COLOR))
483 SKIP(3);
484 }
485 else {
486 if (isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
487 FOR_0_TO_N(3, WRITE(READ()));
488 }
489 else {
490 if (IS_IN_DATA(ATT_COLOR))
491 SKIP(3);
492 }
493 }
494
495 /* set transparency */
496 if (isosurf->att[ATT_TRANSP].changed &&
497 isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
498 if (crnt == 12) {
499 tv = TINTERP(d3, val[ATT_TRANSP]);
500 }
501 else {
502 tv = LINTERP(d, val[ATT_TRANSP][v1], val[ATT_TRANSP][v2]);
503 }
504
505 iso_get_range(isosurf, ATT_TRANSP, &min, &max);
506 c = (min != max) ? 255 - (tv - min) / (max - min) * 255 : 0;
507
508 WRITE(c);
509 if (IS_IN_DATA(ATT_TRANSP))
510 SKIP(1);
511 }
512 else {
513 if (isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
514 WRITE(READ());
515 }
516 else {
517 if (IS_IN_DATA(ATT_TRANSP))
518 SKIP(1);
519 }
520 }
521
522 /* set shin */
523 if (isosurf->att[ATT_SHINE].changed &&
524 isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
525 if (crnt == 12) {
526 tv = TINTERP(d3, val[ATT_SHINE]);
527 }
528 else {
529 tv = LINTERP(d, val[ATT_SHINE][v1], val[ATT_SHINE][v2]);
530 }
531
532 iso_get_range(isosurf, ATT_SHINE, &min, &max);
533 c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
534
535 WRITE(c);
536 if (IS_IN_DATA(ATT_SHINE))
537 SKIP(1);
538 }
539 else {
540 if (isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
541 WRITE(READ());
542 }
543 else {
544 if (IS_IN_DATA(ATT_SHINE))
545 SKIP(1);
546 }
547 }
548
549 /* set emit */
550 if (isosurf->att[ATT_EMIT].changed &&
551 isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
552 if (crnt == 12) {
553 tv = TINTERP(d3, val[ATT_EMIT]);
554 }
555 else {
556 tv = LINTERP(d, val[ATT_EMIT][v1], val[ATT_EMIT][v2]);
557 }
558
559 iso_get_range(isosurf, ATT_EMIT, &min, &max);
560 c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
561
562 WRITE(c);
563 if (IS_IN_DATA(ATT_SHINE))
564 SKIP(1);
565 }
566 else {
567 if (isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
568 WRITE(READ());
569 }
570 else {
571 if (IS_IN_DATA(ATT_EMIT))
572 SKIP(1);
573 }
574 }
575 }
576}
577
578/*!
579 \brief Fill data structure with computed isosurfaces polygons
580
581 \param gvol pointer to geovol struct
582
583 \return 1
584 */
585int gvl_isosurf_calc(geovol *gvol)
586{
587 int x, y, z;
588 int i, a, read;
589 geovol_file *vf;
590 geovol_isosurf *isosurf;
591
592 data_buffer *dbuff;
593 int *need_update, need_update_global;
594
595 dbuff = G_malloc(gvol->n_isosurfs * sizeof(data_buffer));
596 need_update = G_malloc(gvol->n_isosurfs * sizeof(int));
597
598 /* flag - changed any isosurface */
599 need_update_global = 0;
600
601 /* initialize */
602 for (i = 0; i < gvol->n_isosurfs; i++) {
603 isosurf = gvol->isosurf[i];
604
605 /* initialize read/write buffers */
606 dbuff[i].old = NULL;
607 dbuff[i].new = NULL;
608 dbuff[i].ndx_old = 0;
609 dbuff[i].ndx_new = 0;
610 dbuff[i].num_zero = 0;
611
612 need_update[i] = 0;
613 for (a = 1; a < MAX_ATTS; a++) {
614 if (isosurf->att[a].changed) {
615 read = 0;
616 /* changed to map attribute */
617 if (isosurf->att[a].att_src == MAP_ATT) {
618 vf = gvl_file_get_volfile(isosurf->att[a].hfile);
619 read = 1;
620 }
621 /* changed threshold value */
622 if (a == ATT_TOPO) {
623 isosurf->att[a].hfile = gvol->hfile;
624 vf = gvl_file_get_volfile(gvol->hfile);
625 read = 1;
626 }
627 /* initialize reading in selected mode */
628 if (read) {
629 gvl_file_set_mode(vf, 3);
631 }
632
633 /* set update flag - isosurface will be calc */
634 if (read || IS_IN_DATA(a)) {
635 need_update[i] = 1;
636 need_update_global = 1;
637 }
638 }
639 }
640
641 if (need_update[i]) {
642 /* set data buffer */
643 dbuff[i].old = isosurf->data;
644 }
645 }
646
647 /* calculate if only some isosurface changed */
648 if (need_update_global) {
649
650 ResX = gvol->isosurf_x_mod;
651 ResY = gvol->isosurf_y_mod;
652 ResZ = gvol->isosurf_z_mod;
653
654 Cols = gvol->cols / ResX;
655 Rows = gvol->rows / ResY;
656 Depths = gvol->depths / ResZ;
657
658 /* calc isosurface - marching cubes - start */
659
660 for (z = 0; z < Depths - 1; z++) {
661 for (y = 0; y < Rows - 1; y++) {
662 for (x = 0; x < Cols - 1; x++) {
663 for (i = 0; i < gvol->n_isosurfs; i++) {
664 /* recalculate only changed isosurfaces */
665 if (need_update[i]) {
666 iso_calc_cube(gvol->isosurf[i], x, y, z, &dbuff[i]);
667 }
668 }
669 }
670 }
671 }
672 }
673 /* end */
674
675 /* deinitialize */
676 for (i = 0; i < gvol->n_isosurfs; i++) {
677 isosurf = gvol->isosurf[i];
678
679 /* set new isosurface data */
680 if (need_update[i]) {
681 if (dbuff[i].num_zero != 0)
682 gvl_write_char(dbuff[i].ndx_new++, &(dbuff[i].new),
683 dbuff[i].num_zero);
684
685 if (dbuff[i].old == isosurf->data)
686 dbuff[i].old = NULL;
687 G_free(isosurf->data);
688 gvl_align_data(dbuff[i].ndx_new, &(dbuff[i].new));
689 isosurf->data = dbuff[i].new;
690 isosurf->data_desc = 0;
691 }
692
693 for (a = 1; a < MAX_ATTS; a++) {
694 if (isosurf->att[a].changed) {
695 read = 0;
696 /* changed map attribute */
697 if (isosurf->att[a].att_src == MAP_ATT) {
698 vf = gvl_file_get_volfile(isosurf->att[a].hfile);
699 read = 1;
700 }
701 /* changed threshold value */
702 if (a == ATT_TOPO) {
703 isosurf->att[a].hfile = gvol->hfile;
704 vf = gvl_file_get_volfile(gvol->hfile);
705 read = 1;
706 }
707 /* deinitialize reading */
708 if (read) {
710
711 /* set data description */
712 SET_IN_DATA(a);
713 }
714 isosurf->att[a].changed = 0;
715 }
716 else if (isosurf->att[a].att_src == MAP_ATT) {
717 /* set data description */
718 SET_IN_DATA(a);
719 }
720 }
721 }
722
723 /* TODO: G_free() dbuff and need_update ??? */
724
725 return (1);
726}
727
728/*!
729 \brief ADD
730
731 \param pos
732 \param data
733 \param c
734 */
735void gvl_write_char(int pos, unsigned char **data, unsigned char c)
736{
737 /* check to need allocation memory */
738 if ((pos % BUFFER_SIZE) == 0) {
739 *data = G_realloc(*data, sizeof(char) * ((pos / BUFFER_SIZE) + 1) *
741 if (!(*data)) {
742 return;
743 }
744
745 G_debug(3,
746 "gvl_write_char(): reallocate memory for pos : %d to : %lu B",
747 pos, sizeof(char) * ((pos / BUFFER_SIZE) + 1) * BUFFER_SIZE);
748 }
749
750 (*data)[pos] = c;
751}
752
753/*!
754 \brief Read char
755
756 \param pos position index
757 \param data data buffer
758
759 \return char on success
760 \return NULL on failure
761 */
762unsigned char gvl_read_char(int pos, const unsigned char *data)
763{
764 if (!data)
765 return '\0';
766
767 return data[pos];
768}
769
770/*!
771 \brief Append data to buffer
772
773 \param pos position index
774 \param data data buffer
775 */
776void gvl_align_data(int pos, unsigned char **data)
777{
778 unsigned char *p = *data;
779
780 /* realloc memory to fit in data length */
781 p = (unsigned char *)G_realloc(p, sizeof(unsigned char) *
782 pos); /* G_fatal_error */
783 if (!p) {
784 return;
785 }
786
787 G_debug(3, "gvl_align_data(): reallocate memory finally to : %d B", pos);
788
789 if (pos == 0)
790 p = NULL;
791
792 *data = p;
793
794 return;
795}
796
797/************************************************************************/
798/* SLICES */
799
800/************************************************************************/
801
802#define DISTANCE_2(x1, y1, x2, y2) \
803 sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
804
805#define SLICE_MODE_INTERP_NO 0
806#define SLICE_MODE_INTERP_YES 1
807
808/*!
809 \brief Get volume value
810
811 \param gvl pointer to geovol struct
812 \param x,y,z
813
814 \return value
815 */
816float slice_get_value(geovol *gvl, int x, int y, int z)
817{
818 static double d;
819 static geovol_file *vf;
820 static int type;
821 static float value;
822
823 if (x < 0 || y < 0 || z < 0 || (x > gvl->cols - 1) || (y > gvl->rows - 1) ||
824 (z > gvl->depths - 1))
825 return 0.;
826
827 /* get volume file from attribute handle */
828 vf = gvl_file_get_volfile(gvl->hfile);
829 type = gvl_file_get_data_type(vf);
830
831 /* get value from volume file */
832 if (type == VOL_DTYPE_FLOAT) {
833 gvl_file_get_value(vf, x, y, z, &value);
834 }
835 else if (type == VOL_DTYPE_DOUBLE) {
836 gvl_file_get_value(vf, x, y, z, &d);
837 value = (float)d;
838 }
839 else {
840 return 0.;
841 }
842
843 return value;
844}
845
846/*!
847 \brief Calculate slices
848
849 \param gvl pointer to geovol struct
850 \param ndx_slc
851 \param colors
852
853 \return 1
854 */
855int slice_calc(geovol *gvl, int ndx_slc, void *colors)
856{
857 int cols, rows, c, r;
858 int i, j, k, pos, color;
859 int *p_x, *p_y, *p_z;
860 float *p_ex, *p_ey, *p_ez;
861 float value, v[8];
862 float x, y, z, ei, ej, ek, stepx, stepy, stepz;
863 float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
864
865 geovol_slice *slice;
866 geovol_file *vf;
867
868 slice = gvl->slice[ndx_slc];
869
870 /* set mods, pointer to x, y, z step value */
871 if (slice->dir == X) {
872 modx = ResY;
873 mody = ResZ;
874 modz = ResX;
875 p_x = &k;
876 p_y = &i;
877 p_z = &j;
878 p_ex = &ek;
879 p_ey = &ei;
880 p_ez = &ej;
881 }
882 else if (slice->dir == Y) {
883 modx = ResX;
884 mody = ResZ;
885 modz = ResY;
886 p_x = &i;
887 p_y = &k;
888 p_z = &j;
889 p_ex = &ei;
890 p_ey = &ek;
891 p_ez = &ej;
892 }
893 else {
894 modx = ResX;
895 mody = ResY;
896 modz = ResZ;
897 p_x = &i;
898 p_y = &j;
899 p_z = &k;
900 p_ex = &ei;
901 p_ey = &ej;
902 p_ez = &ek;
903 }
904
905 /* distance between slice def. points */
906 distxy = DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
907 distz = fabsf(slice->z2 - slice->z1);
908
909 /* distance between slice def points is zero - nothing to do */
910 if (distxy == 0. || distz == 0.) {
911 return (1);
912 }
913
914 /* start reading volume file */
915 vf = gvl_file_get_volfile(gvl->hfile);
916 gvl_file_set_mode(vf, 3);
918
919 /* set xy resolution */
920 modxy = DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
921 (slice->y2 - slice->y1) / distxy * mody, 0., 0.);
922
923 /* cols/rows of slice */
924 f_cols = distxy / modxy;
925 cols = f_cols > (int)f_cols ? (int)f_cols + 1 : (int)f_cols;
926
927 f_rows = distz / modz;
928 rows = f_rows > (int)f_rows ? (int)f_rows + 1 : (int)f_rows;
929
930 /* set x,y initially to first slice point */
931 x = slice->x1;
932 y = slice->y1;
933
934 /* set x,y step */
935 stepx = (slice->x2 - slice->x1) / f_cols;
936 stepy = (slice->y2 - slice->y1) / f_cols;
937 stepz = (slice->z2 - slice->z1) / f_rows;
938
939 /* set position in slice data */
940 pos = 0;
941
942 /* loop in slice cols */
943 for (c = 0; c < cols + 1; c++) {
944
945 /* convert x, y to integer - index in grid */
946 i = (int)x;
947 j = (int)y;
948
949 /* distance between index and real position */
950 ei = x - (float)i;
951 ej = y - (float)j;
952
953 /* set z to slice z1 point */
954 z = slice->z1;
955
956 /* loop in slice rows */
957 for (r = 0; r < rows + 1; r++) {
958
959 /* distance between index and real position */
960 k = (int)z;
961 ek = z - (float)k;
962
963 /* get interpolated value */
964 if (slice->mode == SLICE_MODE_INTERP_YES) {
965 /* get grid values */
966 v[0] = slice_get_value(gvl, *p_x, *p_y, *p_z);
967 v[1] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z);
968 v[2] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z);
969 v[3] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z);
970
971 v[4] = slice_get_value(gvl, *p_x, *p_y, *p_z + 1);
972 v[5] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z + 1);
973 v[6] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z + 1);
974 v[7] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z + 1);
975
976 /* get interpolated value */
977 value = v[0] * (1. - *p_ex) * (1. - *p_ey) * (1. - *p_ez) +
978 v[1] * (*p_ex) * (1. - *p_ey) * (1. - *p_ez) +
979 v[2] * (1. - *p_ex) * (*p_ey) * (1. - *p_ez) +
980 v[3] * (*p_ex) * (*p_ey) * (1. - *p_ez) +
981 v[4] * (1. - *p_ex) * (1. - *p_ey) * (*p_ez) +
982 v[5] * (*p_ex) * (1. - *p_ey) * (*p_ez) +
983 v[6] * (1. - *p_ex) * (*p_ey) * (*p_ez) +
984 v[7] * (*p_ex) * (*p_ey) * (*p_ez);
985
986 /* no interp value */
987 }
988 else {
989 value = slice_get_value(gvl, *p_x, *p_y, *p_z);
990 }
991
992 /* translate value to color */
993 color = Gvl_get_color_for_value(colors, &value);
994
995 /* write color to slice data */
996 gvl_write_char(pos++, &(slice->data), color & RED_MASK);
997 gvl_write_char(pos++, &(slice->data), (color & GRN_MASK) >> 8);
998 gvl_write_char(pos++, &(slice->data), (color & BLU_MASK) >> 16);
999
1000 /* step in z */
1001 if (r + 1 > f_rows) {
1002 z += stepz * (f_rows - (float)r);
1003 }
1004 else {
1005 z += stepz;
1006 }
1007 }
1008
1009 /* step in x,y */
1010 if (c + 1 > f_cols) {
1011 x += stepx * (f_cols - (float)c);
1012 y += stepy * (f_cols - (float)c);
1013 }
1014 else {
1015 x += stepx;
1016 y += stepy;
1017 }
1018 }
1019
1020 /* end reading volume file */
1022 gvl_align_data(pos, &(slice->data));
1023
1024 return (1);
1025}
1026
1027/*!
1028 \brief Calculate slices for given volume set
1029
1030 \param gvol pointer to geovol struct
1031
1032 \return 1
1033 */
1034int gvl_slices_calc(geovol *gvol)
1035{
1036 int i;
1037 void *colors;
1038
1039 G_debug(5, "gvl_slices_calc(): id=%d", gvol->gvol_id);
1040
1041 /* set current resolution */
1042 ResX = gvol->slice_x_mod;
1043 ResY = gvol->slice_y_mod;
1044 ResZ = gvol->slice_z_mod;
1045
1046 /* set current num of cols, rows, depths */
1047 Cols = gvol->cols / ResX;
1048 Rows = gvol->rows / ResY;
1049 Depths = gvol->depths / ResZ;
1050
1051 /* load colors for geovol file */
1052 Gvl_load_colors_data(&colors, gvl_file_get_name(gvol->hfile));
1053
1054 /* calc changed slices */
1055 for (i = 0; i < gvol->n_slices; i++) {
1056 if (gvol->slice[i]->changed) {
1057 slice_calc(gvol, i, colors);
1058
1059 /* set changed flag */
1060 gvol->slice[i]->changed = 0;
1061 }
1062 }
1063
1064 /* free color */
1065 Gvl_unload_colors_data(colors);
1066
1067 return (1);
1068}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define NULL
Definition ccmath.h:32
CELL_ENTRY cell_table[256]
Definition cell_table.c:3
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double r
int GS_v3norm(float *v1)
Change v1 so that it is a unit vector (2D)
Definition gs_util.c:246
#define RED_MASK
Definition gsd_prim.c:48
#define BLU_MASK
Definition gsd_prim.c:50
#define GRN_MASK
Definition gsd_prim.c:49
int Gvl_unload_colors_data(void *color_data)
Unload color table.
Definition gvl3.c:65
int Gvl_load_colors_data(void **color_data, const char *name)
Load color table.
Definition gvl3.c:34
int Gvl_get_color_for_value(void *color_data, float *value)
Get color for value.
Definition gvl3.c:82
#define DISTANCE_2(x1, y1, x2, y2)
Definition gvl_calc.c:802
int iso_get_cube_values(geovol_isosurf *isosurf, int desc, int x, int y, int z, float *v)
Read values for cube.
Definition gvl_calc.c:228
#define SKIP(n)
Definition gvl_calc.c:60
int mc33_process_cube(int c_ndx, float *v)
ADD.
Definition gvl_calc2.c:307
#define SLICE_MODE_INTERP_YES
Definition gvl_calc.c:806
unsigned char gvl_read_char(int pos, const unsigned char *data)
Read char.
Definition gvl_calc.c:762
double ResZ
Definition gvl_calc.c:80
#define SET_IN_DATA(att)
Definition gvl_calc.c:66
#define TINTERP(d, v)
Definition gvl_calc.c:37
void gvl_align_data(int pos, unsigned char **data)
Append data to buffer.
Definition gvl_calc.c:776
#define BUFFER_SIZE
memory buffer for writing
Definition gvl_calc.c:31
int iso_r_cndx(data_buffer *dbuff)
Read cube index.
Definition gvl_calc.c:126
int slice_calc(geovol *gvl, int ndx_slc, void *colors)
Calculate slices.
Definition gvl_calc.c:855
int gvl_slices_calc(geovol *gvol)
Calculate slices for given volume set.
Definition gvl_calc.c:1034
#define WRITE(c)
writing and reading isosurface data
Definition gvl_calc.c:58
#define FOR_0_TO_N(n, cmd)
Definition gvl_calc.c:47
void iso_w_cndx(int ndx, data_buffer *dbuff)
Write cube index.
Definition gvl_calc.c:91
#define READ()
Definition gvl_calc.c:59
void iso_get_range(geovol_isosurf *isosurf, int desc, double *min, double *max)
Get volume file values range.
Definition gvl_calc.c:212
int Rows
Definition gvl_calc.c:79
#define IS_IN_DATA(att)
check and set data descriptor
Definition gvl_calc.c:65
int Cols
Definition gvl_calc.c:79
void iso_get_cube_grads(geovol_isosurf *isosurf, int x, int y, int z, float(*grad)[3])
Calculate cube grads.
Definition gvl_calc.c:251
void gvl_write_char(int pos, unsigned char **data, unsigned char c)
ADD.
Definition gvl_calc.c:735
int gvl_isosurf_calc(geovol *gvol)
Fill data structure with computed isosurfaces polygons.
Definition gvl_calc.c:585
#define LINTERP(d, a, b)
Definition gvl_calc.c:36
float slice_get_value(geovol *gvl, int x, int y, int z)
Get volume value.
Definition gvl_calc.c:816
double ResX
Definition gvl_calc.c:80
#define FOR_VAR
Definition gvl_calc.c:46
void iso_calc_cube(geovol_isosurf *isosurf, int x, int y, int z, data_buffer *dbuff)
Process cube.
Definition gvl_calc.c:328
double ResY
Definition gvl_calc.c:80
int Depths
Definition gvl_calc.c:79
int iso_get_cube_value(geovol_isosurf *isosurf, int desc, int x, int y, int z, float *v)
Get value from data input.
Definition gvl_calc.c:161
int gvl_file_get_data_type(geovol_file *vf)
Get data type for given handle.
Definition gvl_file.c:202
int gvl_file_end_read(geovol_file *vf)
End read - free buffer memory.
Definition gvl_file.c:1016
void gvl_file_get_min_max(geovol_file *vf, double *min, double *max)
Get minimum and maximum value in volume file.
Definition gvl_file.c:214
int gvl_file_set_mode(geovol_file *vf, IFLAG mode)
Set read mode.
Definition gvl_file.c:1114
int gvl_file_get_value(geovol_file *vf, int x, int y, int z, void *value)
Get value for volume file at x, y, z.
Definition gvl_file.c:1048
int gvl_file_start_read(geovol_file *vf)
Start read - allocate memory buffer a read first data into buffer.
Definition gvl_file.c:965
char * gvl_file_get_name(int id)
Get file name for given handle.
Definition gvl_file.c:165
int gvl_file_is_null_value(geovol_file *vf, void *value)
Check for null value.
Definition gvl_file.c:1085
geovol_file * gvl_file_get_volfile(int id)
Get geovol_file structure for given handle.
Definition gvl_file.c:117
OGSF library -.
#define min(a, b)
#define max(a, b)
#define X(j)
#define x
#define Y(j)