GRASS Programmer's Manual
6.4.2(2012)
|
00001 00019 #include <string.h> 00020 #include <stdlib.h> 00021 00022 #include <grass/gis.h> 00023 #include <grass/gstypes.h> 00024 #include <grass/gsurf.h> 00025 #include <grass/G3d.h> 00026 #include <grass/glocale.h> 00027 00028 #define LUCKY 33 00029 00030 #define MODE_DIRECT 0 00031 #define MODE_SLICE 1 00032 #define MODE_FULL 2 00033 #define MODE_PRELOAD 3 00034 00035 #define MODE_DEFAULT 0 00036 00037 #define STATUS_READY 0 00038 #define STATUS_BUSY 1 00039 00043 typedef struct 00044 { 00045 int num, skip; 00046 int crnt, base; 00047 00048 void *slice[MAX_VOL_SLICES]; 00049 } slice_data; 00050 00051 static geovol_file *Data[MAX_VOL_FILES]; 00052 static geovol_file Df[MAX_VOL_FILES]; /* trying to avoid allocation */ 00053 00054 static int Numfiles = 0; 00055 00056 static int Cur_id = LUCKY; 00057 static int Cur_max; 00058 00059 static int Rows, Cols, Depths; 00060 00061 /* local functions prototypes */ 00062 void *open_g3d_file(const char *, IFLAG *, double *, double *); 00063 int close_g3d_file(void *); 00064 00070 static int init_volfiles(void) 00071 { 00072 int i; 00073 G3D_Region *w3; 00074 00075 for (i = 0; i < MAX_VOL_FILES; i++) { 00076 /* avoiding dynamic allocation */ 00077 Data[i] = &(Df[i]); 00078 } 00079 00080 Cur_max = MAX_VOL_FILES; 00081 00082 /* get window */ 00083 w3 = GVL_get_window(); 00084 00085 /* set cols, rows, depths from window */ 00086 Cols = w3->cols; 00087 Rows = w3->rows; 00088 Depths = w3->depths; 00089 00090 return (1); 00091 } 00092 00098 static int check_num_volfiles(void) 00099 { 00100 if (Numfiles < Cur_max) { 00101 return (0); 00102 } 00103 00104 G_fatal_error(_("Maximum number of datafiles exceeded")); 00105 00106 /* This return statement keeps compilers happy, it is never executed */ 00107 return (0); 00108 } 00109 00118 geovol_file *gvl_file_get_volfile(int id) 00119 { 00120 int i; 00121 00122 for (i = 0; i < Numfiles; i++) { 00123 if (Data[i]->data_id == id) { 00124 return (Data[i]); 00125 } 00126 } 00127 00128 return (NULL); 00129 } 00130 00140 int find_datah(const char *name, IFLAG type, int begin) 00141 { 00142 static int i; 00143 int start; 00144 00145 start = begin ? 0 : i + 1; 00146 00147 for (i = start; i < Numfiles; i++) { 00148 if (!strcmp(Data[i]->file_name, name)) { 00149 if (Data[i]->file_type == type) { 00150 return (Data[i]->data_id); 00151 } 00152 } 00153 } 00154 00155 return (-1); 00156 } 00157 00166 char *gvl_file_get_name(int id) 00167 { 00168 int i; 00169 geovol_file *fvf; 00170 static char retstr[GPATH_MAX]; 00171 00172 for (i = 0; i < Numfiles; i++) { 00173 if (Data[i]->data_id == id) { 00174 fvf = Data[i]; 00175 strcpy(retstr, fvf->file_name); 00176 00177 return (retstr); 00178 } 00179 } 00180 00181 return (NULL); 00182 } 00183 00191 int gvl_file_get_file_type(geovol_file * vf) 00192 { 00193 return (vf->file_type); 00194 } 00195 00203 int gvl_file_get_data_type(geovol_file * vf) 00204 { 00205 return (vf->data_type); 00206 } 00207 00215 void gvl_file_get_min_max(geovol_file * vf, double *min, double *max) 00216 { 00217 *min = vf->min; 00218 *max = vf->max; 00219 } 00220 00233 void *open_volfile(const char *name, IFLAG file_type, IFLAG * data_type, 00234 double *min, double *max) 00235 { 00236 if (file_type == VOL_FTYPE_G3D) { 00237 return open_g3d_file(name, data_type, min, max); 00238 } 00239 00240 return (NULL); 00241 } 00242 00252 int close_volfile(void *map, IFLAG type) 00253 { 00254 if (type == VOL_FTYPE_G3D) { 00255 return close_g3d_file(map); 00256 } 00257 00258 return (-1); 00259 } 00260 00270 int gvl_file_newh(const char *name, IFLAG file_type) 00271 { 00272 geovol_file *new; 00273 static int first = 1; 00274 int i, id; 00275 void *m; 00276 IFLAG data_type; 00277 double min, max; 00278 00279 if (first) { 00280 if (0 > init_volfiles()) { 00281 return (-1); 00282 } 00283 00284 first = 0; 00285 } 00286 00287 if (0 <= (id = find_datah(name, file_type, 1))) { 00288 for (i = 0; i < Numfiles; i++) { 00289 if (Data[i]->data_id == id) { 00290 new = Data[i]; 00291 new->count++; 00292 00293 return (id); 00294 } 00295 } 00296 } 00297 00298 if (0 > check_num_volfiles()) { 00299 return (-1); 00300 } 00301 00302 if (!name) { 00303 return (-1); 00304 } 00305 00306 if ((m = open_volfile(name, file_type, &data_type, &min, &max)) == NULL) { 00307 return (-1); 00308 } 00309 00310 new = Data[Numfiles]; 00311 00312 if (new) { 00313 Numfiles++; 00314 new->data_id = Cur_id++; 00315 00316 new->file_name = G_store(name); 00317 new->file_type = file_type; 00318 new->count = 1; 00319 new->map = m; 00320 new->min = min; 00321 new->max = max; 00322 new->data_type = data_type; 00323 00324 new->status = STATUS_READY; 00325 new->buff = NULL; 00326 00327 new->mode = 255; 00328 gvl_file_set_mode(new, MODE_DEFAULT); 00329 00330 return (new->data_id); 00331 } 00332 00333 return (-1); 00334 } 00335 00343 int free_volfile_buffs(geovol_file * vf) 00344 { 00345 if (vf->mode == MODE_SLICE) { 00346 G_free(vf->buff); 00347 vf->buff = NULL; 00348 } 00349 00350 if (vf->mode == MODE_PRELOAD) { 00351 G_free(vf->buff); 00352 vf->buff = NULL; 00353 } 00354 00355 return (1); 00356 } 00357 00365 int gvl_file_free_datah(int id) 00366 { 00367 int i, j, found = -1; 00368 geovol_file *fvf; 00369 00370 G_debug(5, "gvl_file_free_datah(): id=%d", id); 00371 00372 for (i = 0; i < Numfiles; i++) { 00373 if (Data[i]->data_id == id) { 00374 found = 1; 00375 fvf = Data[i]; 00376 00377 if (fvf->count > 1) { 00378 fvf->count--; 00379 } 00380 else { 00381 close_volfile(fvf->map, fvf->file_type); 00382 free_volfile_buffs(fvf); 00383 00384 G_free(fvf->file_name); 00385 fvf->file_name = NULL; 00386 00387 fvf->data_id = 0; 00388 00389 for (j = i; j < (Numfiles - 1); j++) { 00390 Data[j] = Data[j + 1]; 00391 } 00392 00393 Data[j] = fvf; 00394 00395 --Numfiles; 00396 } 00397 } 00398 } 00399 00400 return (found); 00401 } 00402 00403 /******************************************************************/ 00404 /* reading from G3D raster volume files */ 00405 00406 /******************************************************************/ 00407 00418 void *open_g3d_file(const char *filename, IFLAG * type, double *min, 00419 double *max) 00420 { 00421 const char *mapset; 00422 int itype; 00423 void *map; 00424 00425 /* search for g3d file a return his mapset */ 00426 mapset = G_find_grid3(filename, ""); 00427 if (!mapset) { 00428 G_warning(_("3D raster map <%s> not found"), filename); 00429 return (NULL); 00430 } 00431 00432 /* open g3d file */ 00433 map = 00434 G3d_openCellOld(filename, mapset, G3D_DEFAULT_WINDOW, 00435 G3D_TILE_SAME_AS_FILE, G3D_USE_CACHE_DEFAULT); 00436 if (!map) { 00437 G_warning(_("Unable to open 3D raster map <%s>"), filename); 00438 return (NULL); 00439 } 00440 00441 /* load range into range structure of map */ 00442 if (!G3d_range_load(map)) { 00443 G_warning(_("Unable to read range of 3D raster map <%s>"), filename); 00444 return (NULL); 00445 } 00446 00447 G3d_range_min_max(map, min, max); 00448 00449 /* get file data type */ 00450 itype = G3d_fileTypeMap(map); 00451 if (itype == FCELL_TYPE) 00452 *type = VOL_DTYPE_FLOAT; 00453 if (itype == DCELL_TYPE) 00454 *type = VOL_DTYPE_DOUBLE; 00455 00456 return (map); 00457 } 00458 00467 int close_g3d_file(void *map) 00468 { 00469 /* close opened g3d file */ 00470 if (G3d_closeCell((G3D_Map *) map) != 1) { 00471 G_warning(_("Unable to close 3D raster map <%s>"), 00472 ((G3D_Map *) map)->fileName); 00473 return (-1); 00474 } 00475 00476 return (1); 00477 } 00478 00490 int read_g3d_value(IFLAG type, void *map, int x, int y, int z, void *value) 00491 { 00492 switch (type) { 00493 /* float data type */ 00494 case (VOL_DTYPE_FLOAT): 00495 *((float *)value) = G3d_getFloat(map, x, y, z); 00496 break; 00497 00498 /* double data type */ 00499 case (VOL_DTYPE_DOUBLE): 00500 *((double *)value) = G3d_getDouble(map, x, y, z); 00501 break; 00502 00503 /* unsupported data type */ 00504 default: 00505 return (-1); 00506 } 00507 00508 return (1); 00509 } 00510 00522 int read_g3d_slice(IFLAG type, void *map, int level, void *data) 00523 { 00524 int x, y; 00525 00526 switch (type) { 00527 /* float data type */ 00528 case (VOL_DTYPE_FLOAT): 00529 for (x = 0; x < Cols; x++) { 00530 for (y = 0; y < Rows; y++) { 00531 ((float *)data)[x + y * Cols] = 00532 G3d_getFloat(map, x, y, level); 00533 } 00534 } 00535 00536 break; 00537 00538 /* double data type */ 00539 case (VOL_DTYPE_DOUBLE): 00540 for (x = 0; x < Cols; x++) { 00541 for (y = 0; y < Rows; y++) { 00542 ((double *)data)[x + y * Cols] = 00543 G3d_getDouble(map, x, y, level); 00544 } 00545 } 00546 00547 break; 00548 00549 /* unsupported data type */ 00550 default: 00551 return (-1); 00552 } 00553 00554 return (1); 00555 } 00556 00567 int read_g3d_vol(IFLAG type, void *map, void *data) 00568 { 00569 int x, y, z; 00570 00571 switch (type) { 00572 /* float data type */ 00573 case (VOL_DTYPE_FLOAT): 00574 for (x = 0; x < Cols; x++) { 00575 for (y = 0; y < Rows; y++) { 00576 for (z = 0; z < Depths; z++) { 00577 ((float *)data)[x + y * Cols + z * Rows * Cols] = 00578 G3d_getFloat(map, x, y, z); 00579 } 00580 } 00581 } 00582 00583 break; 00584 00585 /* double data type */ 00586 case (VOL_DTYPE_DOUBLE): 00587 for (x = 0; x < Cols; x++) { 00588 for (y = 0; y < Rows; y++) { 00589 for (z = 0; z < Depths; z++) { 00590 ((double *)data)[x + y * Cols + z * Rows * Cols] = 00591 G3d_getDouble(map, x, y, z); 00592 } 00593 } 00594 } 00595 00596 break; 00597 00598 /* unsupported data type */ 00599 default: 00600 return (-1); 00601 } 00602 00603 return (1); 00604 } 00605 00616 int is_null_g3d_value(IFLAG type, void *value) 00617 { 00618 switch (type) { 00619 /* float data type */ 00620 case (VOL_DTYPE_FLOAT): 00621 return G3d_isNullValueNum(value, FCELL_TYPE); 00622 break; 00623 00624 /* double data type */ 00625 case (VOL_DTYPE_DOUBLE): 00626 return G3d_isNullValueNum(value, DCELL_TYPE); 00627 break; 00628 00629 /* unsupported data type */ 00630 default: 00631 return (-1); 00632 } 00633 00634 return (-1); 00635 } 00636 00637 /******************************************************************/ 00638 /* reading from buffer */ 00639 00640 /******************************************************************/ 00641 00653 int get_buff_value(IFLAG type, void *data, int offset, void *value) 00654 { 00655 switch (type) { 00656 /* float data type */ 00657 case (VOL_DTYPE_FLOAT): 00658 *((float *)value) = ((float *)data)[offset]; 00659 break; 00660 00661 /* double data type */ 00662 case (VOL_DTYPE_DOUBLE): 00663 *((double *)value) = ((double *)data)[offset]; 00664 break; 00665 00666 /* unsupported data type */ 00667 default: 00668 return (-1); 00669 } 00670 00671 return (1); 00672 } 00673 00674 /******************************************************************/ 00675 /* direct mode reading from volume file */ 00676 00677 /******************************************************************/ 00678 00689 int get_direct_value(geovol_file * vf, int x, int y, int z, void *value) 00690 { 00691 switch (vf->file_type) { 00692 /* G3D file type */ 00693 case (VOL_FTYPE_G3D): 00694 if (0 > read_g3d_value(vf->data_type, vf->map, x, y, z, value)) 00695 return (-1); 00696 break; 00697 00698 default: 00699 return (-1); 00700 } 00701 00702 return (1); 00703 } 00704 00705 /******************************************************************/ 00706 /* full mode reading from volume file */ 00707 00708 /******************************************************************/ 00709 00718 int alloc_vol_buff(geovol_file * vf) 00719 { 00720 switch (vf->data_type) { 00721 /* float data type */ 00722 case (VOL_DTYPE_FLOAT): 00723 if ((vf->buff = 00724 (float *)G_malloc(sizeof(float) * Cols * Rows * Depths)) == NULL) 00725 return (-1); 00726 break; 00727 00728 /* double data type */ 00729 case (VOL_DTYPE_DOUBLE): 00730 if ((vf->buff = 00731 (double *)G_malloc(sizeof(double) * Cols * Rows * Depths)) == 00732 NULL) 00733 return (-1); 00734 break; 00735 00736 /* unsupported data type */ 00737 default: 00738 return (-1); 00739 } 00740 00741 return (1); 00742 } 00743 00751 int free_vol_buff(geovol_file * vf) 00752 { 00753 G_free(vf->buff); 00754 00755 return (1); 00756 } 00757 00766 int read_vol(geovol_file * vf) 00767 { 00768 switch (vf->file_type) { 00769 /* G3D file format */ 00770 case (VOL_FTYPE_G3D): 00771 if (0 > read_g3d_vol(vf->data_type, vf->map, vf->buff)) 00772 return (-1); 00773 break; 00774 /* unsupported file format */ 00775 default: 00776 return (-1); 00777 } 00778 00779 return (1); 00780 } 00781 00791 int get_vol_value(geovol_file * vf, int x, int y, int z, void *value) 00792 { 00793 get_buff_value(vf->data_type, vf->buff, z * Rows * Cols + y * Cols + x, 00794 value); 00795 00796 return (1); 00797 } 00798 00799 /******************************************************************/ 00800 /* slice mode reading from volume file */ 00801 00802 /******************************************************************/ 00803 00812 int alloc_slice_buff(geovol_file * vf) 00813 { 00814 int i; 00815 slice_data *sd = (slice_data *) vf->buff; 00816 00817 switch (vf->data_type) { 00818 /* float data type */ 00819 case (VOL_DTYPE_FLOAT): 00820 for (i = 0; i < sd->num; i++) { 00821 if ((sd->slice[i] = 00822 (float *)G_malloc(sizeof(float) * Cols * Rows)) == NULL) 00823 return (-1); 00824 } 00825 break; 00826 00827 /* double data type */ 00828 case (VOL_DTYPE_DOUBLE): 00829 for (i = 0; i < sd->num; i++) { 00830 if ((sd->slice[i] = 00831 (double *)G_malloc(sizeof(double) * Cols * Rows)) == NULL) 00832 return (-1); 00833 } 00834 break; 00835 00836 /* unsupported data type */ 00837 default: 00838 return (-1); 00839 } 00840 00841 return (1); 00842 } 00843 00851 int free_slice_buff(geovol_file * vf) 00852 { 00853 int i; 00854 slice_data *sd = (slice_data *) vf->buff; 00855 00856 for (i = 0; i < sd->num; i++) { 00857 G_free(sd->slice[i]); 00858 } 00859 00860 return (1); 00861 } 00862 00873 int read_slice(geovol_file * vf, int s, int l) 00874 { 00875 /* get slice structure */ 00876 slice_data *sd = (slice_data *) vf->buff; 00877 00878 switch (vf->file_type) { 00879 /* G3D file format */ 00880 case (VOL_FTYPE_G3D): 00881 if (0 > read_g3d_slice(vf->data_type, vf->map, l, sd->slice[s])) 00882 return (-1); 00883 break; 00884 /* unsupported file format */ 00885 default: 00886 return (-1); 00887 } 00888 00889 return (1); 00890 } 00891 00897 void shift_slices(geovol_file * vf) 00898 { 00899 void *tmp; 00900 int i; 00901 00902 slice_data *sd = (slice_data *) vf->buff; 00903 00904 /* change pointers to slices */ 00905 tmp = sd->slice[0]; 00906 for (i = 0; i < sd->num - 1; i++) { 00907 sd->slice[i] = sd->slice[i + 1]; 00908 } 00909 sd->slice[sd->num - 1] = tmp; 00910 00911 /* read new slice data */ 00912 read_slice(vf, sd->num, sd->crnt + 1 + (sd->num - sd->base)); 00913 00914 /* increase current slice value */ 00915 sd->crnt++; 00916 } 00917 00928 int get_slice_value(geovol_file * vf, int x, int y, int z, void *value) 00929 { 00930 slice_data *sd = (slice_data *) vf->buff; 00931 00932 /* value is in loaded slices */ 00933 if ((z >= sd->crnt - (sd->base - 1)) && 00934 (z <= sd->crnt + sd->num - sd->base)) { 00935 00936 get_buff_value(vf->data_type, sd->slice[(z - sd->crnt)], y * Cols + x, 00937 value); 00938 } 00939 00940 /* if value isn't in loaded slices, need read new data slice */ 00941 else if (z == sd->crnt - (sd->base - 1) + 1) { 00942 shift_slices(vf); 00943 get_buff_value(vf->data_type, sd->slice[(z - sd->crnt)], y * Cols + x, 00944 value); 00945 } 00946 00947 /* value out of range */ 00948 else { 00949 return (-1); 00950 } 00951 00952 return (1); 00953 } 00954 00955 /******************************************************************/ 00956 /* reading from volume file */ 00957 00958 /******************************************************************/ 00959 00968 int gvl_file_start_read(geovol_file * vf) 00969 { 00970 int i; 00971 slice_data *sd; 00972 00973 /* check status */ 00974 if (vf->status == STATUS_BUSY) 00975 return (-1); 00976 00977 switch (vf->mode) { 00978 /* read whole volume into memory */ 00979 case (MODE_FULL): 00980 /* allocate memory */ 00981 if (0 > alloc_vol_buff(vf)) 00982 return (-1); 00983 00984 /* read volume */ 00985 read_vol(vf); 00986 break; 00987 00988 /* read first data slices into memory */ 00989 case (MODE_SLICE): 00990 /* allocate slices buffer memory */ 00991 if (0 > alloc_slice_buff(vf)) 00992 return (-1); 00993 00994 /* read volume */ 00995 sd = (slice_data *) vf->buff; 00996 /* set current slice to 0 */ 00997 sd->crnt = 0; 00998 00999 /* read first slices into buffer */ 01000 for (i = 0; i < (sd->num - sd->base + 1); i++) 01001 read_slice(vf, (sd->base - 1) + i, i); 01002 break; 01003 } 01004 01005 /* set status */ 01006 vf->status = STATUS_BUSY; 01007 01008 return (1); 01009 } 01010 01019 int gvl_file_end_read(geovol_file * vf) 01020 { 01021 /* check status */ 01022 if (vf->status == STATUS_READY) 01023 return (-1); 01024 01025 switch (vf->mode) { 01026 case (MODE_FULL): 01027 if (0 > free_vol_buff(vf)) 01028 return (-1); 01029 break; 01030 01031 case (MODE_SLICE): 01032 /* allocate slices buffer memory */ 01033 if (0 > free_slice_buff(vf)) 01034 return (-1); 01035 } 01036 01037 /* set status */ 01038 vf->status = STATUS_READY; 01039 01040 return (1); 01041 } 01042 01051 int gvl_file_get_value(geovol_file * vf, int x, int y, int z, void *value) 01052 { 01053 /* check status */ 01054 if (vf->status != STATUS_BUSY) { 01055 return (-1); 01056 } 01057 01058 switch (vf->mode) { 01059 /* read value direct from g3d file */ 01060 case (MODE_DIRECT): 01061 if (0 > get_direct_value(vf, x, y, z, value)) 01062 return (-1); 01063 break; 01064 01065 case (MODE_SLICE): 01066 if (0 > get_slice_value(vf, x, y, z, value)) 01067 return (-1); 01068 break; 01069 01070 case (MODE_FULL): 01071 case (MODE_PRELOAD): 01072 if (0 > get_vol_value(vf, x, y, z, value)) 01073 return (-1); 01074 break; 01075 } 01076 return (1); 01077 } 01078 01088 int gvl_file_is_null_value(geovol_file * vf, void *value) 01089 { 01090 switch (vf->file_type) { 01091 /* G3D file format */ 01092 case (VOL_FTYPE_G3D): 01093 return is_null_g3d_value(vf->file_type, value); 01094 break; 01095 /* unsupported file format */ 01096 default: 01097 return (-1); 01098 } 01099 01100 return (-1); 01101 } 01102 01103 /******************************************************************/ 01104 /* set parameters for reading volumes */ 01105 01106 /******************************************************************/ 01107 01117 int gvl_file_set_mode(geovol_file * vf, IFLAG mode) 01118 { 01119 slice_data *sd; 01120 01121 if (vf->status == STATUS_BUSY) 01122 return (-1); 01123 01124 if (vf->mode == mode) 01125 return (1); 01126 01127 if (vf->mode == MODE_SLICE) 01128 G_free(vf->buff); 01129 01130 if (vf->mode == MODE_PRELOAD) 01131 G_free(vf->buff); 01132 01133 if (mode == MODE_SLICE) { 01134 if ((vf->buff = (slice_data *) G_malloc(sizeof(slice_data))) == NULL) 01135 return (-1); 01136 01137 sd = (slice_data *) vf->buff; 01138 sd->num = 1; 01139 sd->crnt = 0; 01140 sd->base = 1; 01141 } 01142 01143 if (mode == MODE_PRELOAD) { 01144 /* allocate memory */ 01145 if (0 > alloc_vol_buff(vf)) 01146 return (-1); 01147 01148 /* read volume */ 01149 read_vol(vf); 01150 } 01151 01152 vf->mode = mode; 01153 01154 return (1); 01155 } 01156 01167 int gvl_file_set_slices_param(geovol_file * vf, int n, int b) 01168 { 01169 slice_data *sd; 01170 01171 if (vf->status == STATUS_BUSY) 01172 return (-1); 01173 01174 if (!(vf->mode == MODE_SLICE)) 01175 return (-1); 01176 01177 sd = (slice_data *) vf->buff; 01178 sd->num = n; 01179 sd->base = b; 01180 01181 return (1); 01182 }