GRASS Programmer's Manual  6.4.2(2012)
gvl_calc.c
Go to the documentation of this file.
00001 
00019 #include <math.h>
00020 
00021 #include <grass/gis.h>
00022 #include <grass/gstypes.h>
00023 
00024 #include "rgbpack.h"
00025 #include "mc33_table.h"
00026 
00030 #define BUFFER_SIZE 1000000
00031 
00032 /* USEFUL MACROS */
00033 
00034 /* interp. */
00035 #define LINTERP(d,a,b)  (a + d * (b - a))
00036 #define TINTERP(d,v) ((v[0]*(1.-d[0])*(1.-d[1])*(1.-d[2])) +\
00037                         (v[1]*d[0]*(1.-d[1])*(1.-d[2])) + \
00038                         (v[2]*d[0]*d[1]*(1.-d[2])) + \
00039                         (v[3]*(1.-d[0])*d[1]*(1.-d[2])) + \
00040                         (v[4]*(1.-d[0])*(1.-d[1])*d[2]) + \
00041                         (v[5]*d[0]*(1.-d[1])*d[2]) + \
00042                         (v[6]*d[0]*d[1]*d[2]) + \
00043                         (v[7]*(1.-d[0])*d[1]*d[2]))
00044 
00045 #define FOR_VAR i_for
00046 #define FOR_0_TO_N(n, cmd) { int FOR_VAR; for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) {cmd;} }
00047 
00051 #define WRITE(c) gvl_write_char(dbuff->ndx_new++, &(dbuff->new), c)
00052 #define READ() gvl_read_char(dbuff->ndx_old++, dbuff->old)
00053 #define SKIP(n) dbuff->ndx_old = dbuff->ndx_old + n
00054 
00058 #define IS_IN_DATA(att) ((isosurf->data_desc >> att) & 1)
00059 #define SET_IN_DATA(att) isosurf->data_desc = (isosurf->data_desc | (1 << att))
00060 
00061 typedef struct
00062 {
00063     unsigned char *old;
00064     unsigned char *new;
00065     int ndx_old;
00066     int ndx_new;
00067     int num_zero;
00068 } data_buffer;
00069 
00070 int mc33_process_cube(int c_ndx, float *v);
00071 
00072 /* global variables */
00073 int Rows, Cols, Depths;
00074 double ResX, ResY, ResZ;
00075 
00076 /************************************************************************/
00077 /* ISOSURFACES */
00078 
00085 void iso_w_cndx(int ndx, data_buffer * dbuff)
00086 {
00087     /* cube don't contains polys */
00088     if (ndx == -1) {
00089         if (dbuff->num_zero == 0) {
00090             WRITE(0);
00091             dbuff->num_zero++;
00092         }
00093         else if (dbuff->num_zero == 254) {
00094             WRITE(dbuff->num_zero + 1);
00095             dbuff->num_zero = 0;
00096         }
00097         else {
00098             dbuff->num_zero++;
00099         }
00100     }
00101     else {                      /* isosurface cube */
00102         if (dbuff->num_zero == 0) {
00103             WRITE((ndx / 256) + 1);
00104             WRITE(ndx % 256);
00105         }
00106         else {
00107             WRITE(dbuff->num_zero);
00108             dbuff->num_zero = 0;
00109             WRITE((ndx / 256) + 1);
00110             WRITE(ndx % 256);
00111         }
00112     }
00113 }
00114 
00120 int iso_r_cndx(data_buffer * dbuff)
00121 {
00122     int ndx, ndx2;
00123 
00124     if (dbuff->num_zero != 0) {
00125         dbuff->num_zero--;
00126         ndx = -1;
00127     }
00128     else {
00129         WRITE(ndx = READ());
00130         if (ndx == 0) {
00131             WRITE(dbuff->num_zero = READ());
00132             dbuff->num_zero--;
00133             ndx = -1;
00134         }
00135         else {
00136             WRITE(ndx2 = READ());
00137             ndx = (ndx - 1) * 256 + ndx2;
00138         }
00139     }
00140 
00141     return ndx;
00142 }
00143 
00155 int iso_get_cube_value(geovol_isosurf * isosurf, int desc, int x, int y,
00156                        int z, float *v)
00157 {
00158     double d;
00159     geovol_file *vf;
00160     int type, ret = 1;
00161 
00162     /* get volume file from attribute handle */
00163     vf = gvl_file_get_volfile(isosurf->att[desc].hfile);
00164     type = gvl_file_get_data_type(vf);
00165 
00166     /* get value from volume file */
00167     if (type == VOL_DTYPE_FLOAT) {
00168         gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
00169                            (int)(z * ResZ), v);
00170     }
00171     else if (type == VOL_DTYPE_DOUBLE) {
00172         gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
00173                            (int)(z * ResZ), &d);
00174         *v = (float)d;
00175     }
00176     else {
00177         return 0;
00178     }
00179 
00180     /* null check */
00181     if (gvl_file_is_null_value(vf, v))
00182         ret = 0;
00183 
00184     /* adjust data */
00185     switch (desc) {
00186     case (ATT_TOPO):
00187         *v = (*v) - isosurf->att[desc].constant;
00188         break;
00189     case (ATT_MASK):
00190         if (isosurf->att[desc].constant)
00191             ret = !ret;
00192         break;
00193     }
00194 
00195     return ret;
00196 }
00197 
00206 void iso_get_range(geovol_isosurf * isosurf, int desc, double *min,
00207                    double *max)
00208 {
00209     gvl_file_get_min_max(gvl_file_get_volfile(isosurf->att[desc].hfile), min,
00210                          max);
00211 }
00212 
00223 int iso_get_cube_values(geovol_isosurf * isosurf, int desc, int x, int y,
00224                         int z, float *v)
00225 {
00226     int p, ret = 1;
00227 
00228     for (p = 0; p < 8; ++p) {
00229         if (iso_get_cube_value
00230             (isosurf, desc, x + ((p ^ (p >> 1)) & 1), y + ((p >> 1) & 1),
00231              z + ((p >> 2) & 1), &v[p]) == 0) {
00232             ret = 0;
00233         }
00234     }
00235 
00236     return ret;
00237 }
00238 
00246 void iso_get_cube_grads(geovol_isosurf * isosurf, int x, int y, int z,
00247                         float (*grad)[3])
00248 {
00249     float v[3];
00250     int i, j, k, p;
00251 
00252     for (p = 0; p < 8; ++p) {
00253         i = x + ((p ^ (p >> 1)) & 1);
00254         j = y + ((p >> 1) & 1);
00255         k = z + ((p >> 2) & 1);
00256 
00257         /* x */
00258         if (i == 0) {
00259             iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00260             iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
00261             grad[p][0] = v[2] - v[1];
00262         }
00263         else {
00264             if (i == (Cols - 1)) {
00265                 iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
00266                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00267                 grad[p][0] = v[1] - v[0];
00268             }
00269             else {
00270                 iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
00271                 iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
00272                 grad[p][0] = (v[2] - v[0]) / 2;
00273             }
00274         }
00275 
00276         /* y */
00277         if (j == 0) {
00278             iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00279             iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
00280             grad[p][1] = v[2] - v[1];
00281         }
00282         else {
00283             if (j == (Rows - 1)) {
00284                 iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
00285                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00286                 grad[p][1] = v[1] - v[0];
00287             }
00288             else {
00289                 iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
00290                 iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
00291                 grad[p][1] = (v[2] - v[0]) / 2;
00292             }
00293         }
00294 
00295         /* z */
00296         if (k == 0) {
00297             iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00298             iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
00299             grad[p][2] = v[2] - v[1];
00300         }
00301         else {
00302             if (k == (Depths - 1)) {
00303                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
00304                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
00305                 grad[p][2] = v[1] - v[0];
00306             }
00307             else {
00308                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
00309                 iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
00310                 grad[p][2] = (v[2] - v[0]) / 2;
00311             }
00312         }
00313     }
00314 }
00315 
00323 void iso_calc_cube(geovol_isosurf * isosurf, int x, int y, int z,
00324                    data_buffer * dbuff)
00325 {
00326     int i, c_ndx;
00327     int crnt, v1, v2, c;
00328     float val[MAX_ATTS][8], grad[8][3];
00329     float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
00330     double min, max;
00331 
00332     if (isosurf->att[ATT_TOPO].changed) {
00333         /* read topo values, if there are NULL values then return */
00334         if (!iso_get_cube_values(isosurf, ATT_TOPO, x, y, z, val[ATT_TOPO])) {
00335             iso_w_cndx(-1, dbuff);
00336             return;
00337         }
00338 
00339         /* mask */
00340         if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
00341             if (!iso_get_cube_values
00342                 (isosurf, ATT_MASK, x, y, z, val[ATT_MASK])) {
00343                 iso_w_cndx(-1, dbuff);
00344                 return;
00345             }
00346         }
00347 
00348         /* index to precalculated table */
00349         c_ndx = 0;
00350         for (i = 0; i < 8; i++) {
00351             if (val[ATT_TOPO][i] > 0)
00352                 c_ndx |= 1 << i;
00353         }
00354         c_ndx = mc33_process_cube(c_ndx, val[ATT_TOPO]);
00355 
00356         iso_w_cndx(c_ndx, dbuff);
00357 
00358         if (c_ndx == -1)
00359             return;
00360 
00361         /* calc cube grads */
00362         iso_get_cube_grads(isosurf, x, y, z, grad);
00363 
00364     }
00365     else {
00366         /* read cube index */
00367         if ((c_ndx = iso_r_cndx(dbuff)) == -1)
00368             return;
00369     }
00370 
00371     /* get color values */
00372     if (isosurf->att[ATT_COLOR].changed &&
00373         isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
00374         iso_get_cube_values(isosurf, ATT_COLOR, x, y, z, val[ATT_COLOR]);
00375     }
00376 
00377     /* get transparency values */
00378     if (isosurf->att[ATT_TRANSP].changed &&
00379         isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
00380         iso_get_cube_values(isosurf, ATT_TRANSP, x, y, z, val[ATT_TRANSP]);
00381     }
00382 
00383     /* get shine values */
00384     if (isosurf->att[ATT_SHINE].changed &&
00385         isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
00386         iso_get_cube_values(isosurf, ATT_SHINE, x, y, z, val[ATT_SHINE]);
00387     }
00388 
00389     /* get emit values */
00390     if (isosurf->att[ATT_EMIT].changed &&
00391         isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
00392         iso_get_cube_values(isosurf, ATT_EMIT, x, y, z, val[ATT_EMIT]);
00393     }
00394 
00395     FOR_0_TO_N(3, d_sum[FOR_VAR] = 0.; n_sum[FOR_VAR] = 0.);
00396 
00397     /* loop in edges */
00398     for (i = 0; i < cell_table[c_ndx].nedges; i++) {
00399         /* get edge number */
00400         crnt = cell_table[c_ndx].edges[i];
00401 
00402         /* set topo */
00403         if (isosurf->att[ATT_TOPO].changed) {
00404             /* interior vertex */
00405             if (crnt == 12) {
00406                 FOR_0_TO_N(3,
00407                            WRITE((d3[FOR_VAR] =
00408                                   d_sum[FOR_VAR] /
00409                                   ((float)(cell_table[c_ndx].nedges))) *
00410                                  255));
00411                 GS_v3norm(n_sum);
00412                 FOR_0_TO_N(3,
00413                            WRITE((n_sum[FOR_VAR] /
00414                                   ((float)(cell_table[c_ndx].nedges)) +
00415                                   1.) * 127));
00416                 /* edge vertex */
00417             }
00418             else {
00419                 /* set egdes verts */
00420                 v1 = edge_vert[crnt][0];
00421                 v2 = edge_vert[crnt][1];
00422 
00423                 /* calc intersection point - edge and isosurf */
00424                 d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] -
00425                                          val[ATT_TOPO][v2]);
00426 
00427                 d_sum[edge_vert_pos[crnt][0]] += d;
00428                 d_sum[edge_vert_pos[crnt][1]] += edge_vert_pos[crnt][2];
00429                 d_sum[edge_vert_pos[crnt][3]] += edge_vert_pos[crnt][4];
00430 
00431                 WRITE(d * 255);
00432 
00433                 /* set normal for intersect. point */
00434                 FOR_0_TO_N(3, n[FOR_VAR] =
00435                            LINTERP(d, grad[v1][FOR_VAR], grad[v2][FOR_VAR]));
00436                 GS_v3norm(n);
00437                 FOR_0_TO_N(3, n_sum[FOR_VAR] += n[FOR_VAR]);
00438                 FOR_0_TO_N(3, WRITE((n[FOR_VAR] + 1.) * 127));
00439             }
00440         }
00441         else {
00442             /* read x,y,z of intersection point in cube coords */
00443             if (crnt == 12) {
00444                 WRITE(c = READ());
00445                 d3[0] = ((float)c) / 255.0;
00446                 WRITE(c = READ());
00447                 d3[1] = ((float)c) / 255.0;
00448                 WRITE(c = READ());
00449                 d3[2] = ((float)c) / 255.0;
00450             }
00451             else {
00452                 /* set egdes verts */
00453                 v1 = edge_vert[crnt][0];
00454                 v2 = edge_vert[crnt][1];
00455 
00456                 WRITE(c = READ());
00457                 d = ((float)c) / 255.0;
00458             }
00459 
00460             /* set normals */
00461             FOR_0_TO_N(3, WRITE(READ()));
00462         }
00463 
00464         /* set color */
00465         if (isosurf->att[ATT_COLOR].changed &&
00466             isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
00467             if (crnt == 12) {
00468                 tv = TINTERP(d3, val[ATT_COLOR]);
00469             }
00470             else {
00471                 tv = LINTERP(d, val[ATT_COLOR][v1], val[ATT_COLOR][v2]);
00472             }
00473 
00474             c = Gvl_get_color_for_value(isosurf->att[ATT_COLOR].att_data,
00475                                         &tv);
00476 
00477             WRITE(c & RED_MASK);
00478             WRITE((c & GRN_MASK) >> 8);
00479             WRITE((c & BLU_MASK) >> 16);
00480 
00481             if (IS_IN_DATA(ATT_COLOR))
00482                 SKIP(3);
00483         }
00484         else {
00485             if (isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
00486                 FOR_0_TO_N(3, WRITE(READ()));
00487             }
00488             else {
00489                 if (IS_IN_DATA(ATT_COLOR))
00490                     SKIP(3);
00491             }
00492         }
00493 
00494         /* set transparency */
00495         if (isosurf->att[ATT_TRANSP].changed &&
00496             isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
00497             if (crnt == 12) {
00498                 tv = TINTERP(d3, val[ATT_TRANSP]);
00499             }
00500             else {
00501                 tv = LINTERP(d, val[ATT_TRANSP][v1], val[ATT_TRANSP][v2]);
00502             }
00503 
00504             iso_get_range(isosurf, ATT_TRANSP, &min, &max);
00505             c = (min != max) ? 255 - (tv - min) / (max - min) * 255 : 0;
00506 
00507             WRITE(c);
00508             if (IS_IN_DATA(ATT_TRANSP))
00509                 SKIP(1);
00510         }
00511         else {
00512             if (isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
00513                 WRITE(READ());
00514             }
00515             else {
00516                 if (IS_IN_DATA(ATT_TRANSP))
00517                     SKIP(1);
00518             }
00519         }
00520 
00521         /* set shin */
00522         if (isosurf->att[ATT_SHINE].changed &&
00523             isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
00524             if (crnt == 12) {
00525                 tv = TINTERP(d3, val[ATT_SHINE]);
00526             }
00527             else {
00528                 tv = LINTERP(d, val[ATT_SHINE][v1], val[ATT_SHINE][v2]);
00529             }
00530 
00531             iso_get_range(isosurf, ATT_SHINE, &min, &max);
00532             c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
00533 
00534             WRITE(c);
00535             if (IS_IN_DATA(ATT_SHINE))
00536                 SKIP(1);
00537         }
00538         else {
00539             if (isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
00540                 WRITE(READ());
00541             }
00542             else {
00543                 if (IS_IN_DATA(ATT_SHINE))
00544                     SKIP(1);
00545             }
00546         }
00547 
00548         /* set emit */
00549         if (isosurf->att[ATT_EMIT].changed &&
00550             isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
00551             if (crnt == 12) {
00552                 tv = TINTERP(d3, val[ATT_EMIT]);
00553             }
00554             else {
00555                 tv = LINTERP(d, val[ATT_EMIT][v1], val[ATT_EMIT][v2]);
00556             }
00557 
00558             iso_get_range(isosurf, ATT_EMIT, &min, &max);
00559             c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
00560 
00561             WRITE(c);
00562             if (IS_IN_DATA(ATT_SHINE))
00563                 SKIP(1);
00564         }
00565         else {
00566             if (isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
00567                 WRITE(READ());
00568             }
00569             else {
00570                 if (IS_IN_DATA(ATT_EMIT))
00571                     SKIP(1);
00572             }
00573         }
00574     }
00575 }
00576 
00584 int gvl_isosurf_calc(geovol * gvol)
00585 {
00586     int x, y, z;
00587     int i, a, read;
00588     geovol_file *vf;
00589     geovol_isosurf *isosurf;
00590 
00591     data_buffer *dbuff;
00592     int *need_update, need_update_global;
00593 
00594     dbuff = G_malloc(gvol->n_isosurfs * sizeof(data_buffer));
00595     need_update = G_malloc(gvol->n_isosurfs * sizeof(int));
00596 
00597     /* flag - changed any isosurface */
00598     need_update_global = 0;
00599 
00600     /* initialize */
00601     for (i = 0; i < gvol->n_isosurfs; i++) {
00602         isosurf = gvol->isosurf[i];
00603 
00604         need_update[i] = 0;
00605         for (a = 1; a < MAX_ATTS; a++) {
00606             if (isosurf->att[a].changed) {
00607                 read = 0;
00608                 /* changed to map attribute */
00609                 if (isosurf->att[a].att_src == MAP_ATT) {
00610                     vf = gvl_file_get_volfile(isosurf->att[a].hfile);
00611                     read = 1;
00612                 }
00613                 /* changed threshold value */
00614                 if (a == ATT_TOPO) {
00615                     isosurf->att[a].hfile = gvol->hfile;
00616                     vf = gvl_file_get_volfile(gvol->hfile);
00617                     read = 1;
00618                 }
00619                 /* initialize reading in selected mode */
00620                 if (read) {
00621                     gvl_file_set_mode(vf, 3);
00622                     gvl_file_start_read(vf);
00623                 }
00624 
00625                 /* set update flag - isosurface will be calc */
00626                 if (read || IS_IN_DATA(a)) {
00627                     need_update[i] = 1;
00628                     need_update_global = 1;
00629                 }
00630             }
00631         }
00632 
00633         if (need_update[i]) {
00634             /* initialize read/write buffers */
00635             dbuff[i].old = isosurf->data;
00636             dbuff[i].new = NULL;
00637             dbuff[i].ndx_old = 0;
00638             dbuff[i].ndx_new = 0;
00639             dbuff[i].num_zero = 0;
00640         }
00641     }
00642 
00643     /* calculate if only some isosurface changed */
00644     if (need_update_global) {
00645 
00646         ResX = gvol->isosurf_x_mod;
00647         ResY = gvol->isosurf_y_mod;
00648         ResZ = gvol->isosurf_z_mod;
00649 
00650         Cols = gvol->cols / ResX;
00651         Rows = gvol->rows / ResY;
00652         Depths = gvol->depths / ResZ;
00653 
00654         /* calc isosurface - marching cubes - start */
00655 
00656         for (z = 0; z < Depths - 1; z++) {
00657             for (y = 0; y < Rows - 1; y++) {
00658                 for (x = 0; x < Cols - 1; x++) {
00659                     for (i = 0; i < gvol->n_isosurfs; i++) {
00660                         /* recalculate only changed isosurfaces */
00661                         if (need_update[i]) {
00662                             iso_calc_cube(gvol->isosurf[i], x, y, z,
00663                                           &dbuff[i]);
00664                         }
00665                     }
00666                 }
00667             }
00668         }
00669 
00670     }
00671     /* end */
00672 
00673     /* deinitialize */
00674     for (i = 0; i < gvol->n_isosurfs; i++) {
00675         isosurf = gvol->isosurf[i];
00676 
00677         /* set new isosurface data */
00678         if (need_update[i]) {
00679             if (dbuff[i].num_zero != 0)
00680                 gvl_write_char(dbuff[i].ndx_new++, &(dbuff[i].new),
00681                                dbuff[i].num_zero);
00682 
00683             G_free(isosurf->data);
00684             gvl_align_data(dbuff[i].ndx_new, dbuff[i].new);
00685             isosurf->data = dbuff[i].new;
00686             isosurf->data_desc = 0;
00687         }
00688 
00689         for (a = 1; a < MAX_ATTS; a++) {
00690             if (isosurf->att[a].changed) {
00691                 read = 0;
00692                 /* changed map attribute */
00693                 if (isosurf->att[a].att_src == MAP_ATT) {
00694                     vf = gvl_file_get_volfile(isosurf->att[a].hfile);
00695                     read = 1;
00696                 }
00697                 /* changed threshold value */
00698                 if (a == ATT_TOPO) {
00699                     isosurf->att[a].hfile = gvol->hfile;
00700                     vf = gvl_file_get_volfile(gvol->hfile);
00701                     read = 1;
00702                 }
00703                 /* deinitialize reading */
00704                 if (read) {
00705                     gvl_file_end_read(vf);
00706 
00707                     /* set data description */
00708                     SET_IN_DATA(a);
00709                 }
00710                 isosurf->att[a].changed = 0;
00711             }
00712             else if (isosurf->att[a].att_src == MAP_ATT) {
00713                 /* set data description */
00714                 SET_IN_DATA(a);
00715             }
00716         }
00717     }
00718 
00719     return (1);
00720 }
00721 
00729 void gvl_write_char(int pos, unsigned char **data, unsigned char c)
00730 {
00731     /* check to need allocation memory */
00732     if ((pos % BUFFER_SIZE) == 0) {
00733         *data = (char *)G_realloc(*data,
00734                                   sizeof(char) * ((pos / BUFFER_SIZE) +
00735                                                   1) * BUFFER_SIZE);
00736         if (!(*data)) {
00737             return;
00738         }
00739 
00740         G_debug(3,
00741                 "gvl_write_char(): reallocate memory for pos : %d to : %d B",
00742                 pos, sizeof(char) * ((pos / BUFFER_SIZE) + 1) * BUFFER_SIZE);
00743     }
00744 
00745     (*data)[pos] = c;
00746 }
00747 
00757 unsigned char gvl_read_char(int pos, const unsigned char *data)
00758 {
00759     if (!data)
00760         return '\0';
00761     
00762     return data[pos];
00763 }
00764 
00771 void gvl_align_data(int pos, unsigned char *data)
00772 {
00773     /* realloc memory to fit in data length */
00774     data = (char *)G_realloc(data, sizeof(char) * pos); /* G_fatal_error */
00775     if (!data) {
00776         return;
00777     }
00778 
00779     G_debug(3, "gvl_align_data(): reallocate memory finally to : %d B", pos);
00780 
00781     return;
00782 }
00783 
00784 /************************************************************************/
00785 /* SLICES */
00786 
00787 /************************************************************************/
00788 
00789 #define DISTANCE_2(x1, y1, x2, y2)      sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))
00790 
00791 #define SLICE_MODE_INTERP_NO    0
00792 #define SLICE_MODE_INTERP_YES   1
00793 
00802 float slice_get_value(geovol * gvl, int x, int y, int z)
00803 {
00804     static double d;
00805     static geovol_file *vf;
00806     static int type;
00807     static float value;
00808 
00809     if (x < 0 || y < 0 || z < 0 || (x > gvl->cols - 1) || (y > gvl->rows - 1)
00810         || (z > gvl->depths - 1))
00811         return 0.;
00812 
00813     /* get volume file from attribute handle */
00814     vf = gvl_file_get_volfile(gvl->hfile);
00815     type = gvl_file_get_data_type(vf);
00816 
00817     /* get value from volume file */
00818     if (type == VOL_DTYPE_FLOAT) {
00819         gvl_file_get_value(vf, x, y, z, &value);
00820     }
00821     else if (type == VOL_DTYPE_DOUBLE) {
00822         gvl_file_get_value(vf, x, y, z, &d);
00823         value = (float)d;
00824     }
00825     else {
00826         return 0.;
00827     }
00828 
00829     return value;
00830 }
00831 
00841 int slice_calc(geovol * gvl, int ndx_slc, void *colors)
00842 {
00843     int cols, rows, c, r;
00844     int i, j, k, pos, color;
00845     int *p_x, *p_y, *p_z;
00846     float *p_ex, *p_ey, *p_ez;
00847     float value, v[8];
00848     float x, y, z, ei, ej, ek, stepx, stepy, stepz;
00849     float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
00850 
00851     geovol_slice *slice;
00852     geovol_file *vf;
00853 
00854     slice = gvl->slice[ndx_slc];
00855 
00856     /* set mods, pointer to x, y, z step value */
00857     if (slice->dir == X) {
00858         modx = ResY;
00859         mody = ResZ;
00860         modz = ResX;
00861         p_x = &k;
00862         p_y = &i;
00863         p_z = &j;
00864         p_ex = &ek;
00865         p_ey = &ei;
00866         p_ez = &ej;
00867     }
00868     else if (slice->dir == Y) {
00869         modx = ResX;
00870         mody = ResZ;
00871         modz = ResY;
00872         p_x = &i;
00873         p_y = &k;
00874         p_z = &j;
00875         p_ex = &ei;
00876         p_ey = &ek;
00877         p_ez = &ej;
00878     }
00879     else {
00880         modx = ResX;
00881         mody = ResY;
00882         modz = ResZ;
00883         p_x = &i;
00884         p_y = &j;
00885         p_z = &k;
00886         p_ex = &ei;
00887         p_ey = &ej;
00888         p_ez = &ek;
00889     }
00890 
00891     /* distance between slice def. points */
00892     distxy = DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
00893     distz = fabsf(slice->z2 - slice->z1);
00894 
00895     /* distance between slice def points is zero - nothing to do */
00896     if (distxy == 0. || distz == 0.) {
00897         return (1);
00898     }
00899 
00900     /* start reading volume file */
00901     vf = gvl_file_get_volfile(gvl->hfile);
00902     gvl_file_set_mode(vf, 3);
00903     gvl_file_start_read(vf);
00904 
00905     /* set xy resolution */
00906     modxy =
00907         DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
00908                    (slice->y2 - slice->y1) / distxy * mody, 0., 0.);
00909 
00910     /* cols/rows of slice */
00911     f_cols = distxy / modxy;
00912     cols = f_cols > (int)f_cols ? (int)f_cols + 1 : (int)f_cols;
00913 
00914     f_rows = distz / modz;
00915     rows = f_rows > (int)f_rows ? (int)f_rows + 1 : (int)f_rows;
00916 
00917     /* set x,y intially to first slice point */
00918     x = slice->x1;
00919     y = slice->y1;
00920 
00921     /* set x,y step */
00922     stepx = (slice->x2 - slice->x1) / f_cols;
00923     stepy = (slice->y2 - slice->y1) / f_cols;
00924     stepz = (slice->z2 - slice->z1) / f_rows;
00925 
00926     /* set position in slice data */
00927     pos = 0;
00928 
00929     /* loop in slice cols */
00930     for (c = 0; c < cols + 1; c++) {
00931 
00932         /* convert x, y to integer - index in grid */
00933         i = (int)x;
00934         j = (int)y;
00935 
00936         /* distance between index and real position */
00937         ei = x - (float)i;
00938         ej = y - (float)j;
00939 
00940         /* set z to slice z1 point */
00941         z = slice->z1;
00942 
00943         /* loop in slice rows */
00944         for (r = 0; r < rows + 1; r++) {
00945 
00946             /* distance between index and real position */
00947             k = (int)z;
00948             ek = z - (float)k;
00949 
00950             /* get interpolated value */
00951             if (slice->mode == SLICE_MODE_INTERP_YES) {
00952                 /* get grid values */
00953                 v[0] = slice_get_value(gvl, *p_x, *p_y, *p_z);
00954                 v[1] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z);
00955                 v[2] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z);
00956                 v[3] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z);
00957 
00958                 v[4] = slice_get_value(gvl, *p_x, *p_y, *p_z + 1);
00959                 v[5] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z + 1);
00960                 v[6] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z + 1);
00961                 v[7] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z + 1);
00962 
00963                 /* get interpolated value */
00964                 value = v[0] * (1. - *p_ex) * (1. - *p_ey) * (1. - *p_ez)
00965                     + v[1] * (*p_ex) * (1. - *p_ey) * (1. - *p_ez)
00966                     + v[2] * (1. - *p_ex) * (*p_ey) * (1. - *p_ez)
00967                     + v[3] * (*p_ex) * (*p_ey) * (1. - *p_ez)
00968                     + v[4] * (1. - *p_ex) * (1. - *p_ey) * (*p_ez)
00969                     + v[5] * (*p_ex) * (1. - *p_ey) * (*p_ez)
00970                     + v[6] * (1. - *p_ex) * (*p_ey) * (*p_ez)
00971                     + v[7] * (*p_ex) * (*p_ey) * (*p_ez);
00972 
00973                 /* no interp value */
00974             }
00975             else {
00976                 value = slice_get_value(gvl, *p_x, *p_y, *p_z);
00977             }
00978 
00979             /* translate value to color */
00980             color = Gvl_get_color_for_value(colors, &value);
00981 
00982             /* write color to slice data */
00983             gvl_write_char(pos++, &(slice->data), color & RED_MASK);
00984             gvl_write_char(pos++, &(slice->data), (color & GRN_MASK) >> 8);
00985             gvl_write_char(pos++, &(slice->data), (color & BLU_MASK) >> 16);
00986 
00987             /* step in z */
00988             if (r + 1 > f_rows) {
00989                 z += stepz * (f_rows - (float)r);
00990             }
00991             else {
00992                 z += stepz;
00993             }
00994         }
00995 
00996         /* step in x,y */
00997         if (c + 1 > f_cols) {
00998             x += stepx * (f_cols - (float)c);
00999             y += stepy * (f_cols - (float)c);
01000         }
01001         else {
01002             x += stepx;
01003             y += stepy;
01004         }
01005     }
01006 
01007     /* end reading volume file */
01008     gvl_file_end_read(vf);
01009     gvl_align_data(pos, slice->data);
01010 
01011     return (1);
01012 }
01013 
01021 int gvl_slices_calc(geovol *gvol)
01022 {
01023     int i;
01024     void *colors;
01025 
01026     G_debug(5, "gvl_slices_calc(): id=%d", gvol->gvol_id);
01027     
01028     /* set current resolution */
01029     ResX = gvol->slice_x_mod;
01030     ResY = gvol->slice_y_mod;
01031     ResZ = gvol->slice_z_mod;
01032 
01033     /* set current num of cols, rows, depths */
01034     Cols = gvol->cols / ResX;
01035     Rows = gvol->rows / ResY;
01036     Depths = gvol->depths / ResZ;
01037 
01038     /* load colors for geovol file */
01039     Gvl_load_colors_data(&colors, gvl_file_get_name(gvol->hfile));
01040 
01041     /* calc changed slices */
01042     for (i = 0; i < gvol->n_slices; i++) {
01043         if (gvol->slice[i]->changed) {
01044             slice_calc(gvol, i, colors);
01045 
01046             /* set changed flag */
01047             gvol->slice[i]->changed = 0;
01048         }
01049     }
01050 
01051     /* free color */
01052     Gvl_unload_colors_data(colors);
01053 
01054     return (1);
01055 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines