GRASS Programmer's Manual
6.4.2(2012)
|
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 }