GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 
61 typedef 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 
70 int mc33_process_cube(int c_ndx, float *v);
71 
72 /* global variables */
74 double ResX, ResY, ResZ;
75 
76 /************************************************************************/
77 /* ISOSURFACES */
78 
79 /*!
80  \brief Write cube index
81 
82  \param ndx
83  \param dbuff
84  */
85 void 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  */
120 int 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  */
155 int 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  */
206 void 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  */
223 int 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  */
246 void 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  */
323 void 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 
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  */
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) {
711  gvl_file_end_read(vf);
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  */
737 void 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  */
765 unsigned 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  */
779 void 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  */
818 float 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  */
857 int 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 */
1024  gvl_file_end_read(vf);
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  */
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 }
int Rows
Definition: gvl_calc.c:73
int slice_x_mod
Definition: ogsf.h:455
#define G_malloc(n)
Definition: defs/gis.h:112
OGSF library -.
#define BLU_MASK
Definition: ogsf.h:198
int data_desc
Definition: ogsf.h:420
void gvl_file_get_min_max(geovol_file *, double *, double *)
Get minimum and maximum value in volume file.
Definition: gvl_file.c:214
int rows
Definition: ogsf.h:440
#define ATT_SHINE
Definition: ogsf.h:77
#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_z_mod
Definition: ogsf.h:455
#define RED_MASK
Definition: ogsf.h:196
#define IS_IN_DATA(att)
check and set data descriptor
Definition: gvl_calc.c:58
float constant
Definition: ogsf.h:409
IFLAG att_src
Definition: ogsf.h:405
int Gvl_unload_colors_data(void *)
Unload color table.
Definition: gvl3.c:65
int depths
Definition: ogsf.h:440
#define LINTERP(d, a, b)
Definition: gvl_calc.c:35
#define min(x, y)
Definition: draw2.c:31
int gvl_file_get_value(geovol_file *, int, int, int, void *)
Get value for volume file at x, y, z.
Definition: gvl_file.c:1050
#define ATT_TOPO
Definition: ogsf.h:73
int gvl_file_is_null_value(geovol_file *, void *)
Check for null value.
Definition: gvl_file.c:1087
double ResX
Definition: gvl_calc.c:74
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
geovol_isosurf * isosurf[MAX_ISOSURFS]
Definition: ogsf.h:449
#define FOR_0_TO_N(n, cmd)
Definition: gvl_calc.c:46
#define NULL
Definition: ccmath.h:32
int mc33_process_cube(int c_ndx, float *v)
ADD.
Definition: gvl_calc2.c:307
#define x
#define max(x, y)
Definition: draw2.c:32
#define VOL_DTYPE_FLOAT
Definition: ogsf.h:131
int gvl_isosurf_calc(geovol *gvol)
Fill data structure with computed isosurfaces polygons.
Definition: gvl_calc.c:585
double ResY
Definition: gvl_calc.c:74
int GS_v3norm(float *)
Change v1 so that it is a unit vector (2D)
Definition: gs_util.c:246
int Gvl_load_colors_data(void **, const char *)
Load color table.
Definition: gvl3.c:34
#define TINTERP(d, v)
Definition: gvl_calc.c:36
float z2
Definition: ogsf.h:427
int isosurf_y_mod
Definition: ogsf.h:450
#define SET_IN_DATA(att)
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:206
int n_isosurfs
Definition: ogsf.h:448
int edges[12]
Definition: viz.h:76
int gvl_file_get_data_type(geovol_file *)
Get data type for given handle.
Definition: gvl_file.c:202
float y2
Definition: ogsf.h:427
int cols
Definition: ogsf.h:440
void gvl_write_char(int pos, unsigned char **data, unsigned char c)
ADD.
Definition: gvl_calc.c:737
int n_slices
Definition: ogsf.h:453
int gvol_id
Definition: ogsf.h:436
int isosurf_x_mod
Definition: ogsf.h:450
geovol_isosurf_att att[MAX_ATTS]
Definition: ogsf.h:418
int nedges
Definition: viz.h:75
unsigned char gvl_read_char(int pos, const unsigned char *data)
Read char.
Definition: gvl_calc.c:765
#define SKIP(n)
Definition: gvl_calc.c:53
#define MAP_ATT
Definition: ogsf.h:83
char * gvl_file_get_name(int)
Get file name for given handle.
Definition: gvl_file.c:165
geovol_slice * slice[MAX_SLICES]
Definition: ogsf.h:454
#define ATT_COLOR
Definition: ogsf.h:74
int dir
Definition: ogsf.h:426
void iso_calc_cube(geovol_isosurf *isosurf, int x, int y, int z, data_buffer *dbuff)
Process cube.
Definition: gvl_calc.c:323
int Gvl_get_color_for_value(void *, float *)
Get color for value.
Definition: gvl3.c:82
void iso_w_cndx(int ndx, data_buffer *dbuff)
Write cube index.
Definition: gvl_calc.c:85
unsigned char * data
Definition: ogsf.h:421
#define SLICE_MODE_INTERP_YES
Definition: gvl_calc.c:808
#define WRITE(c)
writing and reading isosurface data
Definition: gvl_calc.c:51
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
Definition: ogsf.h:434
#define Y
Definition: ogsf.h:138
int slice_y_mod
Definition: ogsf.h:455
int hfile
Definition: ogsf.h:439
void * att_data
Definition: ogsf.h:411
void gvl_align_data(int pos, unsigned char **data)
Append data to buffer.
Definition: gvl_calc.c:779
#define ATT_TRANSP
Definition: ogsf.h:76
unsigned char * data
Definition: ogsf.h:428
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
geovol_file * gvl_file_get_volfile(int)
Get geovol_file structure for given handle.
Definition: gvl_file.c:117
float slice_get_value(geovol *gvl, int x, int y, int z)
Get volume value.
Definition: gvl_calc.c:818
#define DISTANCE_2(x1, y1, x2, y2)
Definition: gvl_calc.c:805
int isosurf_z_mod
Definition: ogsf.h:450
#define G_realloc(p, n)
Definition: defs/gis.h:114
int gvl_file_end_read(geovol_file *)
End read - free buffer memory.
Definition: gvl_file.c:1018
#define X
Definition: ogsf.h:137
int gvl_file_set_mode(geovol_file *, IFLAG)
Set read mode.
Definition: gvl_file.c:1116
#define MAX_ATTS
Definition: ogsf.h:43
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
int mode
Definition: ogsf.h:431
float z1
Definition: ogsf.h:427
#define ATT_MASK
Definition: ogsf.h:75
float x2
Definition: ogsf.h:427
#define VOL_DTYPE_DOUBLE
Definition: ogsf.h:132
float x1
Definition: ogsf.h:427
#define ATT_EMIT
Definition: ogsf.h:78
int Cols
Definition: gvl_calc.c:73
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
CELL_ENTRY cell_table[256]
Definition: cell_table.c:3
int G_debug(int, const char *,...) __attribute__((format(printf
int gvl_file_start_read(geovol_file *)
Start read - allocate memory buffer a read first data into buffer.
Definition: gvl_file.c:967
#define READ()
Definition: gvl_calc.c:52
double r
Definition: r_raster.c:39
float y1
Definition: ogsf.h:427
#define FOR_VAR
Definition: gvl_calc.c:45
int changed
Definition: ogsf.h:429
#define GRN_MASK
Definition: ogsf.h:197
double ResZ
Definition: gvl_calc.c:74
int Depths
Definition: gvl_calc.c:73