21 #include <grass/gis.h>
22 #include <grass/ogsf.h>
30 #define BUFFER_SIZE 1000000
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]))
46 #define FOR_0_TO_N(n, cmd) { int FOR_VAR; for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) {cmd;} }
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
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))
89 if (dbuff->num_zero == 0) {
93 else if (dbuff->num_zero == 254) {
94 WRITE(dbuff->num_zero + 1);
102 if (dbuff->num_zero == 0) {
103 WRITE((ndx / 256) + 1);
107 WRITE(dbuff->num_zero);
109 WRITE((ndx / 256) + 1);
124 if (dbuff->num_zero != 0) {
137 ndx = (ndx - 1) * 256 + ndx2;
167 if (type == VOL_DTYPE_FLOAT) {
171 else if (type == VOL_DTYPE_DOUBLE) {
173 (
int)(z *
ResZ), &d);
187 *v = (*v) - isosurf->att[desc].constant;
190 if (isosurf->att[desc].constant)
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) {
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);
261 grad[p][0] = v[2] - v[1];
264 if (i == (
Cols - 1)) {
267 grad[p][0] = v[1] - v[0];
272 grad[p][0] = (v[2] - v[0]) / 2;
280 grad[p][1] = v[2] - v[1];
283 if (j == (
Rows - 1)) {
286 grad[p][1] = v[1] - v[0];
291 grad[p][1] = (v[2] - v[0]) / 2;
299 grad[p][2] = v[2] - v[1];
305 grad[p][2] = v[1] - v[0];
310 grad[p][2] = (v[2] - v[0]) / 2;
328 float val[MAX_ATTS][8], grad[8][3];
329 float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
332 if (isosurf->att[ATT_TOPO].changed) {
340 if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
342 (isosurf, ATT_MASK,
x, y, z, val[ATT_MASK])) {
350 for (i = 0; i < 8; i++) {
351 if (val[ATT_TOPO][i] > 0)
372 if (isosurf->att[ATT_COLOR].changed &&
373 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
378 if (isosurf->att[ATT_TRANSP].changed &&
379 isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
384 if (isosurf->att[ATT_SHINE].changed &&
385 isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
390 if (isosurf->att[ATT_EMIT].changed &&
391 isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
404 if (isosurf->att[ATT_TOPO].changed) {
421 v1 = edge_vert[crnt][0];
422 v2 = edge_vert[crnt][1];
425 d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] -
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];
446 d3[0] = ((float)c) / 255.0;
448 d3[1] = ((float)c) / 255.0;
450 d3[2] = ((float)c) / 255.0;
454 v1 = edge_vert[crnt][0];
455 v2 = edge_vert[crnt][1];
458 d = ((float)c) / 255.0;
466 if (isosurf->att[ATT_COLOR].changed &&
467 isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
469 tv =
TINTERP(d3, val[ATT_COLOR]);
472 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]) {
678 for (i = 0; i < gvol->n_isosurfs; i++) {
679 isosurf = gvol->isosurf[i];
682 if (need_update[i]) {
683 if (dbuff[i].num_zero != 0)
687 if (dbuff[i].old == isosurf->data)
691 isosurf->data = dbuff[i].new;
692 isosurf->data_desc = 0;
695 for (a = 1; a < MAX_ATTS; a++) {
696 if (isosurf->att[a].changed) {
699 if (isosurf->att[a].att_src == MAP_ATT) {
705 isosurf->att[a].hfile = gvol->hfile;
716 isosurf->att[a].changed = 0;
718 else if (isosurf->att[a].att_src == MAP_ATT) {
741 *data = (
char *)G_realloc(*data,
749 "gvl_write_char(): reallocate memory for pos : %d to : %lu B",
781 unsigned char *p = *data;
785 p = (
unsigned char *)G_realloc(p,
sizeof(
unsigned char) * pos);
790 G_debug(3,
"gvl_align_data(): reallocate memory finally to : %d B", pos);
805 #define DISTANCE_2(x1, y1, x2, y2) sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))
807 #define SLICE_MODE_INTERP_NO 0
808 #define SLICE_MODE_INTERP_YES 1
821 static geovol_file *vf;
825 if (
x < 0 || y < 0 || z < 0 || (
x > gvl->cols - 1) || (y > gvl->rows - 1)
826 || (z > gvl->depths - 1))
834 if (type == VOL_DTYPE_FLOAT) {
837 else if (type == VOL_DTYPE_DOUBLE) {
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;
864 float x, y, z, ei, ej, ek, stepx, stepy, stepz;
865 float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
870 slice = gvl->slice[ndx_slc];
873 if (slice->dir ==
X) {
884 else if (slice->dir ==
Y) {
908 distxy =
DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
909 distz = fabsf(slice->z2 - slice->z1);
912 if (distxy == 0. || distz == 0.) {
923 DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
924 (slice->y2 - slice->y1) / distxy * mody, 0., 0.);
927 f_cols = distxy / modxy;
928 cols = f_cols > (int)f_cols ? (
int)f_cols + 1 : (int)f_cols;
930 f_rows = distz / modz;
931 rows = f_rows > (int)f_rows ? (
int)f_rows + 1 : (int)f_rows;
938 stepx = (slice->x2 - slice->x1) / f_cols;
939 stepy = (slice->y2 - slice->y1) / f_cols;
940 stepz = (slice->z2 - slice->z1) / f_rows;
946 for (c = 0; c < cols + 1; c++) {
960 for (
r = 0;
r < rows + 1;
r++) {
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);
1004 if (
r + 1 > f_rows) {
1005 z += stepz * (f_rows - (float)
r);
1013 if (c + 1 > f_cols) {
1014 x += stepx * (f_cols - (float)c);
1015 y += stepy * (f_cols - (float)c);
1042 G_debug(5,
"gvl_slices_calc(): id=%d", gvol->gvol_id);
1045 ResX = gvol->slice_x_mod;
1046 ResY = gvol->slice_y_mod;
1047 ResZ = gvol->slice_z_mod;
1058 for (i = 0; i < gvol->n_slices; i++) {
1059 if (gvol->slice[i]->changed) {
1063 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.
geovol_file * gvl_file_get_volfile(int id)
Get geovol_file structure for given handle.
char * gvl_file_get_name(int id)
Get file name for given handle.
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.
int gvl_file_is_null_value(geovol_file *vf, void *value)
Check for null value.