23#include <grass/ogsf.h>
31#define BUFFER_SIZE 1000000
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]))
47#define FOR_0_TO_N(n, cmd) \
50 for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) { \
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
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))
95 if (dbuff->num_zero == 0) {
99 else if (dbuff->num_zero == 254) {
100 WRITE(dbuff->num_zero + 1);
108 if (dbuff->num_zero == 0) {
109 WRITE((ndx / 256) + 1);
113 WRITE(dbuff->num_zero);
115 WRITE((ndx / 256) + 1);
130 if (dbuff->num_zero != 0) {
143 ndx = (ndx - 1) * 256 + ndx2;
173 if (type == VOL_DTYPE_FLOAT) {
177 else if (type == VOL_DTYPE_DOUBLE) {
179 (
int)(z *
ResZ), &d);
193 *v = (*v) - isosurf->att[desc].constant;
196 if (isosurf->att[desc].constant)
233 for (p = 0; p < 8; ++p) {
235 y + ((p >> 1) & 1), z + ((p >> 2) & 1),
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);
266 grad[p][0] = v[2] - v[1];
269 if (i == (
Cols - 1)) {
272 grad[p][0] = v[1] - v[0];
277 grad[p][0] = (v[2] - v[0]) / 2;
285 grad[p][1] = v[2] - v[1];
288 if (j == (
Rows - 1)) {
291 grad[p][1] = v[1] - v[0];
296 grad[p][1] = (v[2] - v[0]) / 2;
304 grad[p][2] = v[2] - v[1];
310 grad[p][2] = v[1] - v[0];
315 grad[p][2] = (v[2] - v[0]) / 2;
333 float val[MAX_ATTS][8], grad[8][3];
334 float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
337 if (isosurf->att[ATT_TOPO].changed) {
345 if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
355 for (i = 0; i < 8; i++) {
356 if (val[ATT_TOPO][i] > 0)
376 if (isosurf->att[ATT_COLOR].changed &&
377 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
382 if (isosurf->att[ATT_TRANSP].changed &&
383 isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
388 if (isosurf->att[ATT_SHINE].changed &&
389 isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
394 if (isosurf->att[ATT_EMIT].changed &&
395 isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
402 for (i = 0; i <
cell_table[c_ndx].nedges; i++) {
407 if (isosurf->att[ATT_TOPO].changed) {
423 v1 = edge_vert[crnt][0];
424 v2 = edge_vert[crnt][1];
427 d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] - val[ATT_TOPO][v2]);
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];
447 d3[0] = ((float)c) / 255.0;
449 d3[1] = ((float)c) / 255.0;
451 d3[2] = ((float)c) / 255.0;
455 v1 = edge_vert[crnt][0];
456 v2 = edge_vert[crnt][1];
459 d = ((float)c) / 255.0;
467 if (isosurf->att[ATT_COLOR].changed &&
468 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
470 tv =
TINTERP(d3, val[ATT_COLOR]);
473 tv =
LINTERP(d, val[ATT_COLOR][v1], val[ATT_COLOR][v2]);
486 if (isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
496 if (isosurf->att[ATT_TRANSP].changed &&
497 isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
499 tv =
TINTERP(d3, val[ATT_TRANSP]);
502 tv =
LINTERP(d, val[ATT_TRANSP][v1], val[ATT_TRANSP][v2]);
513 if (isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
523 if (isosurf->att[ATT_SHINE].changed &&
524 isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
526 tv =
TINTERP(d3, val[ATT_SHINE]);
529 tv =
LINTERP(d, val[ATT_SHINE][v1], val[ATT_SHINE][v2]);
540 if (isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
550 if (isosurf->att[ATT_EMIT].changed &&
551 isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
553 tv =
TINTERP(d3, val[ATT_EMIT]);
556 tv =
LINTERP(d, val[ATT_EMIT][v1], val[ATT_EMIT][v2]);
567 if (isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
590 geovol_isosurf *isosurf;
593 int *need_update, need_update_global;
595 dbuff = G_malloc(gvol->n_isosurfs *
sizeof(data_buffer));
596 need_update = G_malloc(gvol->n_isosurfs *
sizeof(
int));
599 need_update_global = 0;
602 for (i = 0; i < gvol->n_isosurfs; i++) {
603 isosurf = gvol->isosurf[i];
608 dbuff[i].ndx_old = 0;
609 dbuff[i].ndx_new = 0;
610 dbuff[i].num_zero = 0;
613 for (a = 1; a < MAX_ATTS; a++) {
614 if (isosurf->att[a].changed) {
617 if (isosurf->att[a].att_src == MAP_ATT) {
623 isosurf->att[a].hfile = gvol->hfile;
636 need_update_global = 1;
641 if (need_update[i]) {
643 dbuff[i].old = isosurf->data;
648 if (need_update_global) {
650 ResX = gvol->isosurf_x_mod;
651 ResY = gvol->isosurf_y_mod;
652 ResZ = gvol->isosurf_z_mod;
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++) {
665 if (need_update[i]) {
676 for (i = 0; i < gvol->n_isosurfs; i++) {
677 isosurf = gvol->isosurf[i];
680 if (need_update[i]) {
681 if (dbuff[i].num_zero != 0)
685 if (dbuff[i].old == isosurf->data)
689 isosurf->data = dbuff[i].new;
690 isosurf->data_desc = 0;
693 for (a = 1; a < MAX_ATTS; a++) {
694 if (isosurf->att[a].changed) {
697 if (isosurf->att[a].att_src == MAP_ATT) {
703 isosurf->att[a].hfile = gvol->hfile;
714 isosurf->att[a].changed = 0;
716 else if (isosurf->att[a].att_src == MAP_ATT) {
739 *data = G_realloc(*data,
sizeof(
char) * ((pos /
BUFFER_SIZE) + 1) *
746 "gvl_write_char(): reallocate memory for pos : %d to : %lu B",
778 unsigned char *p = *data;
781 p = (
unsigned char *)G_realloc(p,
sizeof(
unsigned char) *
787 G_debug(3,
"gvl_align_data(): reallocate memory finally to : %d B", pos);
802#define DISTANCE_2(x1, y1, x2, y2) \
803 sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
805#define SLICE_MODE_INTERP_NO 0
806#define SLICE_MODE_INTERP_YES 1
819 static geovol_file *vf;
823 if (
x < 0 || y < 0 || z < 0 || (
x > gvl->cols - 1) || (y > gvl->rows - 1) ||
824 (z > gvl->depths - 1))
832 if (type == VOL_DTYPE_FLOAT) {
835 else if (type == VOL_DTYPE_DOUBLE) {
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;
862 float x, y, z, ei, ej, ek, stepx, stepy, stepz;
863 float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
868 slice = gvl->slice[ndx_slc];
871 if (slice->dir ==
X) {
882 else if (slice->dir ==
Y) {
906 distxy =
DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
907 distz = fabsf(slice->z2 - slice->z1);
910 if (distxy == 0. || distz == 0.) {
920 modxy =
DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
921 (slice->y2 - slice->y1) / distxy * mody, 0., 0.);
924 f_cols = distxy / modxy;
925 cols = f_cols > (int)f_cols ? (
int)f_cols + 1 : (int)f_cols;
927 f_rows = distz / modz;
928 rows = f_rows > (int)f_rows ? (
int)f_rows + 1 : (int)f_rows;
935 stepx = (slice->x2 - slice->x1) / f_cols;
936 stepy = (slice->y2 - slice->y1) / f_cols;
937 stepz = (slice->z2 - slice->z1) / f_rows;
943 for (c = 0; c < cols + 1; c++) {
957 for (
r = 0;
r < rows + 1;
r++) {
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);
1001 if (
r + 1 > f_rows) {
1002 z += stepz * (f_rows - (float)
r);
1010 if (c + 1 > f_cols) {
1011 x += stepx * (f_cols - (float)c);
1012 y += stepy * (f_cols - (float)c);
1039 G_debug(5,
"gvl_slices_calc(): id=%d", gvol->gvol_id);
1042 ResX = gvol->slice_x_mod;
1043 ResY = gvol->slice_y_mod;
1044 ResZ = gvol->slice_z_mod;
1055 for (i = 0; i < gvol->n_slices; i++) {
1056 if (gvol->slice[i]->changed) {
1060 gvol->slice[i]->changed = 0;
void G_free(void *buf)
Free allocated memory.
CELL_ENTRY cell_table[256]
int G_debug(int level, const char *msg,...)
Print debugging message.
int GS_v3norm(float *v1)
Change v1 so that it is a unit vector (2D)
int Gvl_unload_colors_data(void *color_data)
Unload color table.
int Gvl_load_colors_data(void **color_data, const char *name)
Load color table.
int Gvl_get_color_for_value(void *color_data, float *value)
Get color for value.
#define DISTANCE_2(x1, y1, x2, y2)
int iso_get_cube_values(geovol_isosurf *isosurf, int desc, int x, int y, int z, float *v)
Read values for cube.
int mc33_process_cube(int c_ndx, float *v)
ADD.
#define SLICE_MODE_INTERP_YES
unsigned char gvl_read_char(int pos, const unsigned char *data)
Read char.
void gvl_align_data(int pos, unsigned char **data)
Append data to buffer.
#define BUFFER_SIZE
memory buffer for writing
int iso_r_cndx(data_buffer *dbuff)
Read cube index.
int slice_calc(geovol *gvl, int ndx_slc, void *colors)
Calculate slices.
int gvl_slices_calc(geovol *gvol)
Calculate slices for given volume set.
#define WRITE(c)
writing and reading isosurface data
#define FOR_0_TO_N(n, cmd)
void iso_w_cndx(int ndx, data_buffer *dbuff)
Write cube index.
void iso_get_range(geovol_isosurf *isosurf, int desc, double *min, double *max)
Get volume file values range.
#define IS_IN_DATA(att)
check and set data descriptor
void iso_get_cube_grads(geovol_isosurf *isosurf, int x, int y, int z, float(*grad)[3])
Calculate cube grads.
void gvl_write_char(int pos, unsigned char **data, unsigned char c)
ADD.
int gvl_isosurf_calc(geovol *gvol)
Fill data structure with computed isosurfaces polygons.
float slice_get_value(geovol *gvl, int x, int y, int z)
Get volume value.
void iso_calc_cube(geovol_isosurf *isosurf, int x, int y, int z, data_buffer *dbuff)
Process cube.
int iso_get_cube_value(geovol_isosurf *isosurf, int desc, int x, int y, int z, float *v)
Get value from data input.
int gvl_file_get_data_type(geovol_file *vf)
Get data type for given handle.
int gvl_file_end_read(geovol_file *vf)
End read - free buffer memory.
void gvl_file_get_min_max(geovol_file *vf, double *min, double *max)
Get minimum and maximum value in volume file.
int gvl_file_set_mode(geovol_file *vf, IFLAG mode)
Set read mode.
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.
int gvl_file_start_read(geovol_file *vf)
Start read - allocate memory buffer a read first data into buffer.
char * gvl_file_get_name(int id)
Get file name for given handle.
int gvl_file_is_null_value(geovol_file *vf, void *value)
Check for null value.
geovol_file * gvl_file_get_volfile(int id)
Get geovol_file structure for given handle.