diff --git a/src/gmt_api.c b/src/gmt_api.c index 9790f9792c9..46fc1a44c4e 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -2828,6 +2828,107 @@ bool gmtlib_data_is_geographic (struct GMTAPI_CTRL *API, const char *file) { return (geo); } +/*! . */ +GMT_LOCAL unsigned int gmtapi_expand_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *new_wesn, unsigned int *orig_pad, double *orig_wesn) { + unsigned int tmp_pad[4] = {0, 0, 0, 0}, delta[4] = {0, 0, 0, 0}, k = 0; + /* When using subset with memory grids we cannot actually cut the grid but instead + * must temporarily change the pad to match the desired inner region wesn. This means + * the pads will change and can be quite large. */ + struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h); + + gmt_M_memcpy (tmp_pad, h->pad, 4, unsigned int); /* Initialize new pad to the original pad */ + /* First determine which (and how many, k) of the 4 new boundaries are inside the original region and update the padding: */ + if (new_wesn[XLO] > h->wesn[XLO]) k++, tmp_pad[XLO] += urint ((new_wesn[XLO] - h->wesn[XLO]) * HH->r_inc[GMT_X]); + if (new_wesn[XHI] < h->wesn[XHI]) k++, tmp_pad[XHI] += urint ((h->wesn[XHI] - new_wesn[XHI]) * HH->r_inc[GMT_X]); + if (new_wesn[YLO] > h->wesn[YLO]) k++, tmp_pad[YLO] += urint ((new_wesn[YLO] - h->wesn[YLO]) * HH->r_inc[GMT_Y]); + if (new_wesn[YHI] < h->wesn[YHI]) k++, tmp_pad[YHI] += urint ((h->wesn[YHI] - new_wesn[YHI]) * HH->r_inc[GMT_Y]); + if (k) { /* Yes, pad will change since region is different for k of the 4 sides */ + for (k = 0; k < 4; k++) delta[k] = tmp_pad[k] - h->pad[k]; /* Columns with data being passed as padding */ + gmt_M_memcpy (orig_pad, h->pad, 4, unsigned int); /* Place the original grid pad in the provided array */ + gmt_M_memcpy (orig_wesn, h->wesn, 4, double); /* Place the original grid wesn in the provided array */ + gmt_M_memcpy (h->pad, tmp_pad, 4, unsigned int); /* Place the new pad in the grid header */ + gmt_M_memcpy (h->wesn, new_wesn, 4, double); /* Place the new wesn in the grid header */ + gmt_set_grddim (GMT, h); /* This recomputes n_columns|n_rows. */ + GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtapi_expand_headerpad: %d pad sides changed. Now %u/%u/%u/%u\n", + k, h->pad[XLO], h->pad[XHI], h->pad[YLO], h->pad[YHI]); + for (k = 0; k < 4; k++) { /* If pad now contains data then change the BC to reflect this */ + if (delta[k] >= orig_pad[k]) HH->BC[k] = GMT_BC_IS_DATA; + } + } + else + GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtapi_expand_headerpad: No pad adjustment needed\n"); + return k; +} + +/*! . */ +GMT_LOCAL void gmtapi_update_grid_minmax (struct GMT_CTRL *GMT, struct GMT_GRID *G) { + /* Update grid header z_min/z_max to reflect the range within the subset */ + uint64_t ij; + unsigned int row, col; + struct GMT_GRID_HEADER *h = G->header; + struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h); + + h->z_min = +DBL_MAX; h->z_max = -DBL_MAX; /* Reset the min/max before we search */ + HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */ + gmt_M_grd_loop (GMT, G, row, col, ij) { + if (gmt_M_is_fnan (G->data[ij])) + HH->has_NaNs = GMT_GRID_HAS_NANS; + else { + h->z_min = MIN (h->z_min, G->data[ij]); + h->z_max = MAX (h->z_max, G->data[ij]); + } + } +} + +/*! . */ +GMT_LOCAL void gmtapi_update_cube_minmax (struct GMT_CTRL *GMT, struct GMT_CUBE *U) { + /* Update cube header z_min/z_max to reflect the range within the subset */ + uint64_t ij, node, here = 0; + unsigned int k, row, col; + struct GMT_GRID_HEADER *h = U->header; + struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h); + + h->z_min = +DBL_MAX; h->z_max = -DBL_MAX; /* Reset the min/max before we search */ + HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */ + for (k = 0; k < h->n_bands; k++) { + gmt_M_grd_loop (GMT, U, row, col, ij) { + node = ij + here; + if (gmt_M_is_fnan (U->data[node])) + HH->has_NaNs = GMT_GRID_HAS_NANS; + else { + h->z_min = MIN (h->z_min, U->data[node]); + h->z_max = MAX (h->z_max, U->data[node]); + } + } + here += h->size; + } +} + +/*! . */ +GMT_LOCAL void gmtapi_contract_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, unsigned int *orig_pad, double *orig_wesn) { + /* When using subset with memory grids we must reset the pad back to the original setting when done */ + if (h == NULL) return; /* Nothing for us to work with */ + gmt_M_memcpy (h->pad, orig_pad, 4, unsigned int); /* Place the original pad in the grid header */ + gmt_M_memcpy (h->wesn, orig_wesn, 4, double); /* Place the orig_pad wesn in the grid header */ + gmt_set_grddim (GMT, h); /* This recomputes n_columns|n_rows. */ + GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtapi_contract_headerpad: Pad and wesn reset to original values\n"); +} + +/*! . */ +GMT_LOCAL void gmtapi_contract_pad (struct GMT_CTRL *GMT, void *object, int family, unsigned int *orig_pad, double *orig_wesn) { + /* When using subset with memory grids we must reset the pad back to the original setting when done */ + struct GMT_GRID_HEADER *h = NULL; + if (family == GMT_IS_GRID) { + struct GMT_GRID *G = gmtapi_get_grid_data (object); + if (G) h = G->header; + } + else if (family == GMT_IS_IMAGE) { + struct GMT_IMAGE *I = gmtapi_get_image_data (object); + if (I) h = I->header; + } + gmtapi_contract_headerpad (GMT, h, orig_pad, orig_wesn); +} + /*! . */ GMT_LOCAL int gmtapi_get_object (struct GMTAPI_CTRL *API, int sfamily, void *ptr) { /* Returns the ID of the first object whose resource pointer matches ptr. @@ -2879,15 +2980,17 @@ GMT_LOCAL void * gmtapi_pass_object (struct GMTAPI_CTRL *API, struct GMTAPI_DATA } } else if (object->region) { /* Possibly adjust the pad so inner region matches requested wesn */ - /* NOTE: This assumes the memory cannot be adjusted. Probably should distinaguish between GMT_IS_REFERENCE and GMT_IS_DUPLICATE + /* NOTE: This assumes the memory cannot be adjusted. Probably should distinguish between GMT_IS_REFERENCE and GMT_IS_DUPLICATE * and offer different behavior. As it is we assume read-only grids */ if (object->reset_pad) { /* First undo any prior sub-region used with this memory grid */ - gmtlib_contract_headerpad (API->GMT, G->header, object->orig_pad, object->orig_wesn); + gmtapi_contract_headerpad (API->GMT, G->header, object->orig_pad, object->orig_wesn); object->reset_pad = HH->reset_pad = 0; } /* Then apply the new pad adjustment. Basically we cannot mess with the data so we change what constitute the pad */ - if (gmtlib_expand_headerpad (API->GMT, G->header, object->wesn, object->orig_pad, object->orig_wesn)) + if (gmtapi_expand_headerpad (API->GMT, G->header, object->wesn, object->orig_pad, object->orig_wesn)) { object->reset_pad = HH->reset_pad = 1; + gmtapi_update_grid_minmax (API->GMT, G); /* Update z-range */ + } } } if (mode & GMT_CONTAINER_ONLY) break; /* No grid yet */ @@ -2915,11 +3018,11 @@ GMT_LOCAL void * gmtapi_pass_object (struct GMTAPI_CTRL *API, struct GMTAPI_DATA /* NOTE: This assumes the memory cannot be adjusted. Probably should distinaguish between GMT_IS_REFERENCE and GMT_IS_DUPLICATE * and offer different behavior. As it is we assume read-only images */ if (object->reset_pad) { /* First undo a prior sub-region used with this memory grid */ - gmtlib_contract_headerpad (API->GMT, I->header, object->orig_pad, object->orig_wesn); + gmtapi_contract_headerpad (API->GMT, I->header, object->orig_pad, object->orig_wesn); object->reset_pad = HH->reset_pad = 0; } /* Then apply the new pad adjustment. Basically we cannot mess with the data so we change what constitute the pad */ - if (gmtlib_expand_headerpad (API->GMT, I->header, object->wesn, object->orig_pad, object->orig_wesn)) + if (gmtapi_expand_headerpad (API->GMT, I->header, object->wesn, object->orig_pad, object->orig_wesn)) object->reset_pad = HH->reset_pad = 1; } } @@ -4615,10 +4718,10 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob if (done && S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ HH = gmt_get_H_hidden (I_obj->header); if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory image */ - gmtlib_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = HH->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) S_obj->reset_pad = HH->reset_pad = 1; } GMT_Report (API, GMT_MSG_DEBUG, "gmtapi_import_image: Return from GMT_IS_REFERENCE\n"); @@ -4659,10 +4762,10 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob via = true; if (S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory image */ - gmtlib_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = HH->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) S_obj->reset_pad = HH->reset_pad = 1; } break; @@ -4701,10 +4804,10 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob via = true; if (S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory image */ - gmtlib_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, I_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = HH->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, I_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) S_obj->reset_pad = HH->reset_pad = 1; } break; @@ -5118,11 +5221,13 @@ GMT_LOCAL struct GMT_GRID *gmtapi_import_grid (struct GMTAPI_CTRL *API, int obje if (done && S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ HH = gmt_get_H_hidden (G_obj->header); if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory grid */ - gmtlib_contract_headerpad (GMT, G_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, G_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = HH->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, G_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, G_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) { S_obj->reset_pad = HH->reset_pad = 1; + gmtapi_update_grid_minmax (API->GMT, G_obj); /* Update z-range */ + } } GMT_Report (API, GMT_MSG_DEBUG, "gmtapi_import_grid: Return from GMT_IS_REFERENCE\n"); break; @@ -5318,11 +5423,13 @@ GMT_LOCAL struct GMT_GRID *gmtapi_import_grid (struct GMTAPI_CTRL *API, int obje } else if (S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory grid */ - gmtlib_contract_headerpad (GMT, G_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, G_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, G_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, G_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) { S_obj->reset_pad = 1; + gmtapi_update_grid_minmax (GMT, G_obj); /* Update z-range */ + } } break; @@ -6126,11 +6233,13 @@ GMT_LOCAL struct GMT_CUBE * gmtapi_import_cube (struct GMTAPI_CTRL *API, int obj } else if (S_obj->region) { /* Possibly adjust the pad so inner region matches wesn */ if (S_obj->reset_pad) { /* First undo a prior sub-region used with this memory cube */ - gmtlib_contract_headerpad (GMT, U_obj->header, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_headerpad (GMT, U_obj->header, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = 0; } - if (gmtlib_expand_headerpad (GMT, U_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) + if (gmtapi_expand_headerpad (GMT, U_obj->header, S_obj->wesn, S_obj->orig_pad, S_obj->orig_wesn)) { S_obj->reset_pad = 1; + gmtapi_update_cube_minmax (API->GMT, U_obj); /* Update z-range */ + } } break; @@ -7711,7 +7820,7 @@ void gmtlib_garbage_collection (struct GMTAPI_CTRL *API, int level) { if (!(level == GMT_NOTSET || S_obj->alloc_level == u_level)) { /* Not the right module level (or not end of session yet) */ if (S_obj->reset_pad && S_obj->no_longer_owner == false) { /* Temporarily changed pad to access a sub-region of a memory grid - now reset this if still the owner */ address = S_obj->resource; /* Try to get the data object */ - gmtlib_contract_pad (API->GMT, address, S_obj->actual_family, S_obj->orig_pad, S_obj->orig_wesn); + gmtapi_contract_pad (API->GMT, address, S_obj->actual_family, S_obj->orig_pad, S_obj->orig_wesn); S_obj->reset_pad = 0; } i++; continue; diff --git a/src/gmt_grdio.c b/src/gmt_grdio.c index d147a22ea91..a7d0c81885c 100644 --- a/src/gmt_grdio.c +++ b/src/gmt_grdio.c @@ -94,9 +94,6 @@ struct GRD_PAD { /* Local structure */ /* Local functions */ -GMT_LOCAL inline struct GMT_GRID * gmtgrdio_get_grid_data (struct GMT_GRID *ptr) {return (ptr);} -GMT_LOCAL inline struct GMT_IMAGE * gmtgrdio_get_image_data (struct GMT_IMAGE *ptr) {return (ptr);} - /*! gmt_M_grd_get_size computes grid size including the padding, and doubles it if complex values */ GMT_LOCAL size_t gmtgrdio_grd_get_size (struct GMT_GRID_HEADER *h) { return ((((h->complex_mode & GMT_GRID_IS_COMPLEX_MASK) > 0) + 1ULL) * h->mx * h->my); @@ -3232,60 +3229,6 @@ bool gmtlib_found_url_for_gdal (char *fname) { return false; } -unsigned int gmtlib_expand_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *new_wesn, unsigned int *orig_pad, double *orig_wesn) { - unsigned int tmp_pad[4] = {0, 0, 0, 0}, delta[4] = {0, 0, 0, 0}, k = 0; - /* When using subset with memory grids we cannot actually cut the grid but instead - * must temporarily change the pad to match the desired inner region wesn. This means - * the pads will change and can be quite large. */ - struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h); - - gmt_M_memcpy (tmp_pad, h->pad, 4, unsigned int); /* Initialize new pad to the original pad */ - /* First determine which (and how many, k) of the 4 new boundaries are inside the original region and update the padding: */ - if (new_wesn[XLO] > h->wesn[XLO]) k++, tmp_pad[XLO] += urint ((new_wesn[XLO] - h->wesn[XLO]) * HH->r_inc[GMT_X]); - if (new_wesn[XHI] < h->wesn[XHI]) k++, tmp_pad[XHI] += urint ((h->wesn[XHI] - new_wesn[XHI]) * HH->r_inc[GMT_X]); - if (new_wesn[YLO] > h->wesn[YLO]) k++, tmp_pad[YLO] += urint ((new_wesn[YLO] - h->wesn[YLO]) * HH->r_inc[GMT_Y]); - if (new_wesn[YHI] < h->wesn[YHI]) k++, tmp_pad[YHI] += urint ((h->wesn[YHI] - new_wesn[YHI]) * HH->r_inc[GMT_Y]); - if (k) { /* Yes, pad will change since region is different for k of the 4 sides */ - for (k = 0; k < 4; k++) delta[k] = tmp_pad[k] - h->pad[k]; /* Columns with data being passed as padding */ - gmt_M_memcpy (orig_pad, h->pad, 4, unsigned int); /* Place the original grid pad in the provided array */ - gmt_M_memcpy (orig_wesn, h->wesn, 4, double); /* Place the original grid wesn in the provided array */ - gmt_M_memcpy (h->pad, tmp_pad, 4, unsigned int); /* Place the new pad in the grid header */ - gmt_M_memcpy (h->wesn, new_wesn, 4, double); /* Place the new wesn in the grid header */ - gmt_set_grddim (GMT, h); /* This recomputes n_columns|n_rows. */ - GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtlib_expand_headerpad: %d pad sides changed. Now %u/%u/%u/%u\n", - k, h->pad[XLO], h->pad[XHI], h->pad[YLO], h->pad[YHI]); - for (k = 0; k < 4; k++) { /* If pad now contains data then change the BC to reflect this */ - if (delta[k] >= orig_pad[k]) HH->BC[k] = GMT_BC_IS_DATA; - } - } - else - GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtlib_expand_headerpad: No pad adjustment needed\n"); - return k; -} - -void gmtlib_contract_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, unsigned int *orig_pad, double *orig_wesn) { - /* When using subset with memory grids we must reset the pad back to the original setting when done */ - if (h == NULL) return; /* Nothing for us to work with */ - gmt_M_memcpy (h->pad, orig_pad, 4, unsigned int); /* Place the original pad in the grid header */ - gmt_M_memcpy (h->wesn, orig_wesn, 4, double); /* Place the orig_pad wesn in the grid header */ - gmt_set_grddim (GMT, h); /* This recomputes n_columns|n_rows. */ - GMT_Report (GMT->parent, GMT_MSG_DEBUG, "gmtlib_contract_headerpad: Pad and wesn reset to original values\n"); -} - -void gmtlib_contract_pad (struct GMT_CTRL *GMT, void *object, int family, unsigned int *orig_pad, double *orig_wesn) { - /* When using subset with memory grids we must reset the pad back to the original setting when done */ - struct GMT_GRID_HEADER *h = NULL; - if (family == GMT_IS_GRID) { - struct GMT_GRID *G = gmtgrdio_get_grid_data (object); - if (G) h = G->header; - } - else if (family == GMT_IS_IMAGE) { - struct GMT_IMAGE *I = gmtgrdio_get_image_data (object); - if (I) h = I->header; - } - gmtlib_contract_headerpad (GMT, h, orig_pad, orig_wesn); -} - GMT_LOCAL int gmtgrdio_get_extension_period (char *file) { int i, pos_ext = 0; for (i = (int)strlen(file) - 1; i > 0; i--) { diff --git a/src/gmt_internals.h b/src/gmt_internals.h index 756a1f1659b..b21189a4c9e 100644 --- a/src/gmt_internals.h +++ b/src/gmt_internals.h @@ -278,8 +278,6 @@ EXTERN_MSC unsigned int gmtlib_get_arc (struct GMT_CTRL *GMT, double x0, double EXTERN_MSC struct GMT_PALETTE * gmtlib_read_cpt (struct GMT_CTRL *GMT, void *source, unsigned int source_type, unsigned int cpt_flags); EXTERN_MSC int gmtlib_alloc_univector (struct GMT_CTRL *GMT, union GMT_UNIVECTOR *u, unsigned int type, uint64_t n_rows); EXTERN_MSC unsigned int gmtlib_get_arc (struct GMT_CTRL *GMT, double x0, double y0, double r, double dir1, double dir2, double **x, double **y); -EXTERN_MSC unsigned int gmtlib_expand_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double *new_wesn, unsigned int *orig_pad, double *orig_wesn); -EXTERN_MSC void gmtlib_contract_headerpad (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, unsigned int *orig_pad, double *orig_wesn); EXTERN_MSC void gmtlib_contract_pad (struct GMT_CTRL *GMT, void *object, int family, unsigned int *orig_pad, double *orig_wesn); EXTERN_MSC uint64_t gmtlib_glob_list (struct GMT_CTRL *GMT, const char *pattern, char ***list); EXTERN_MSC void gmtlib_change_out_dataset (struct GMT_CTRL *GMT, struct GMT_DATASET *D);