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