GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2653d7b383
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 
68 typedef 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 
76 int mc33_process_cube(int c_ndx, float *v);
77 
78 /* global variables */
80 double ResX, ResY, ResZ;
81 
82 /************************************************************************/
83 /* ISOSURFACES */
84 
85 /*!
86  \brief Write cube index
87 
88  \param ndx
89  \param dbuff
90  */
91 void 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  */
126 int 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  */
161 int 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  */
212 void iso_get_range(geovol_isosurf *isosurf, int desc, double *min, double *max)
213 {
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  */
228 int 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  */
251 void 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  */
328 void 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  */
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) {
709  gvl_file_end_read(vf);
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  */
735 void 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) *
740  BUFFER_SIZE);
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  */
762 unsigned 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  */
776 void 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  */
816 float 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  */
855 int 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 */
1021  gvl_file_end_read(vf);
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  */
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 }
#define NULL
Definition: ccmath.h:32
CELL_ENTRY cell_table[256]
Definition: cell_table.c:3
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_realloc(p, n)
Definition: defs/gis.h:96
#define G_malloc(n)
Definition: defs/gis.h:94
int G_debug(int, const char *,...) __attribute__((format(printf
int GS_v3norm(float *)
Change v1 so that it is a unit vector (2D)
Definition: gs_util.c:246
int gvl_file_end_read(geovol_file *)
End read - free buffer memory.
Definition: gvl_file.c:1016
int gvl_file_start_read(geovol_file *)
Start read - allocate memory buffer a read first data into buffer.
Definition: gvl_file.c:965
char * gvl_file_get_name(int)
Get file name for given handle.
Definition: gvl_file.c:165
int gvl_file_get_data_type(geovol_file *)
Get data type for given handle.
Definition: gvl_file.c:202
int Gvl_get_color_for_value(void *, float *)
Get color for value.
Definition: gvl3.c:82
int gvl_file_set_mode(geovol_file *, IFLAG)
Set read mode.
Definition: gvl_file.c:1114
int Gvl_load_colors_data(void **, const char *)
Load color table.
Definition: gvl3.c:34
void gvl_file_get_min_max(geovol_file *, double *, double *)
Get minimum and maximum value in volume file.
Definition: gvl_file.c:214
int gvl_file_is_null_value(geovol_file *, void *)
Check for null value.
Definition: gvl_file.c:1085
int gvl_file_get_value(geovol_file *, int, int, int, void *)
Get value for volume file at x, y, z.
Definition: gvl_file.c:1048
geovol_file * gvl_file_get_volfile(int)
Get geovol_file structure for given handle.
Definition: gvl_file.c:117
int Gvl_unload_colors_data(void *)
Unload color table.
Definition: gvl3.c:65
#define min(x, y)
Definition: draw2.c:29
#define max(x, y)
Definition: draw2.c:30
#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
OGSF library -.
#define ATT_MASK
Definition: ogsf.h:77
#define VOL_DTYPE_FLOAT
Definition: ogsf.h:134
#define RED_MASK
Definition: ogsf.h:200
#define X
Definition: ogsf.h:140
#define MAX_ATTS
Definition: ogsf.h:45
#define ATT_TOPO
Definition: ogsf.h:75
#define BLU_MASK
Definition: ogsf.h:202
#define ATT_COLOR
Definition: ogsf.h:76
#define ATT_EMIT
Definition: ogsf.h:80
#define GRN_MASK
Definition: ogsf.h:201
#define ATT_SHINE
Definition: ogsf.h:79
#define Y
Definition: ogsf.h:141
#define MAP_ATT
Definition: ogsf.h:85
#define ATT_TRANSP
Definition: ogsf.h:78
#define VOL_DTYPE_DOUBLE
Definition: ogsf.h:135
double r
Definition: r_raster.c:39
int nedges
Definition: viz.h:68
int edges[12]
Definition: viz.h:69
Definition: ogsf.h:426
int n_slices
Definition: ogsf.h:444
int hfile
Definition: ogsf.h:430
int depths
Definition: ogsf.h:431
geovol_slice * slice[MAX_SLICES]
Definition: ogsf.h:445
int slice_z_mod
Definition: ogsf.h:446
int isosurf_z_mod
Definition: ogsf.h:441
int slice_y_mod
Definition: ogsf.h:446
int isosurf_y_mod
Definition: ogsf.h:441
int cols
Definition: ogsf.h:431
int isosurf_x_mod
Definition: ogsf.h:441
int gvol_id
Definition: ogsf.h:427
int slice_x_mod
Definition: ogsf.h:446
geovol_isosurf * isosurf[MAX_ISOSURFS]
Definition: ogsf.h:440
int n_isosurfs
Definition: ogsf.h:439
int rows
Definition: ogsf.h:431
void * att_data
Definition: ogsf.h:405
float constant
Definition: ogsf.h:403
IFLAG att_src
Definition: ogsf.h:399
unsigned char * data
Definition: ogsf.h:414
geovol_isosurf_att att[MAX_ATTS]
Definition: ogsf.h:411
int data_desc
Definition: ogsf.h:413
unsigned char * data
Definition: ogsf.h:420
float z1
Definition: ogsf.h:419
float x1
Definition: ogsf.h:419
float z2
Definition: ogsf.h:419
float y1
Definition: ogsf.h:419
int changed
Definition: ogsf.h:421
int dir
Definition: ogsf.h:418
int mode
Definition: ogsf.h:423
float x2
Definition: ogsf.h:419
float y2
Definition: ogsf.h:419
#define x