/*-------------------------------------------------------------------- * * Copyright (c) 1991-2020 by the GMT Team (https://www.generic-mapping-tools.org/team.html) * See LICENSE.TXT file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation; version 3 or any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * Contact info: www.generic-mapping-tools.org *--------------------------------------------------------------------*/ /* * * G M T _ S U P P O R T . C * *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * GMT_support.c contains code used by most GMT programs * * Author: Paul Wessel * Date: 1-JAN-2010 * Version: 5 * * Modules in this file: * * gmtlib_akima Akima's 1-D spline * gmt_BC_init Initialize BCs for a grid or image * gmt_grd_BC_set Set two rows of padding according to bound cond for grid * gmtlib_image_BC_set Set two rows of padding according to bound cond for image * gmtsupport_check_rgb Check rgb for valid range * gmtsupport_cmyk_to_rgb Corvert CMYK to RGB * gmtsupport_comp_double_asc Used when sorting doubles into ascending order [checks for NaN] * gmtsupport_comp_float_asc Used when sorting floats into ascending order [checks for NaN] * gmtsupport_comp_int_asc Used when sorting ints into ascending order * gmt_contours Subroutine for contouring * gmtlib_cspline Natural cubic 1-D spline solver * gmtsupport_csplint Natural cubic 1-D spline evaluator * gmt_delaunay Performs a Delaunay triangulation * gmtlib_get_annot_label Construct degree/minute label * gmt_get_index Return color table entry for given z * gmt_get_fill_from_z Return fill type for given z * gmt_get_format Find # of decimals and create format string * gmt_get_rgb_from_z Return rgb for given z * gmt_get_plot_array Allocate memory for plotting arrays * gmt_getfill Decipher and check fill argument * gmt_getinc Decipher and check increment argument * gmt_getpen Decipher and check pen argument * gmt_getrgb Decipher and check color argument * gmt_hsv_to_rgb Convert HSV to RGB * gmt_init_fill Initialize fill attributes * gmt_init_pen Initialize pen attributes * gmt_illuminate Add illumination effects to rgb * gmt_intpol 1-D interpolation * gmt_lab_to_rgb Corvert CIELAB LAB to RGB * gmt_lab_to_xyz Convert CIELAB LAB to XYZ * gmt_non_zero_winding Finds if a point is inside/outside a polygon * gmt_putpen Encode pen argument into textstring * gmtlib_read_cpt Read color palette file * gmtsupport_rgb_to_cmyk Convert RGB to CMYK * gmt_rgb_to_hsv Convert RGB to HSV * gmt_rgb_to_lab Convert RGB to CMYK * gmt_rgb_to_xyz Convert RGB to CIELAB XYZ * gmt_sample_cpt Resamples the current CPT based on new z-array * gmt_invert_cpt Flips the current CPT upside down * gmtsupport_smooth_contour Use Akima's spline to smooth contour * GMT_shift_refpoint Adjust reference point based on size and justification of plotted item * gmt_sprintf_float Make formatted string from float, while checking for %-apostrophe * gmtsupport_trace_contour Function that trace the contours in gmt_contours * gmtsupport_polar_adjust Adjust label justification for polar projection * gmt_xyz_to_rgb Convert CIELAB XYZ to RGB * gmt_xyz_to_lab Convert CIELAB XYZ to LAB */ /*! * \file gmt_support.c * \brief GMT_support.c contains code used by most GMT programs. */ #include "gmt_dev.h" #include "gmt_internals.h" #include #ifndef WIN32 #include #else #include #include #endif /*! . */ enum GMT_profmode { GMT_GOT_AZIM = 1, GMT_GOT_ORIENT = 2, GMT_GOT_LENGTH = 4, GMT_GOT_NP = 8, GMT_GOT_INC = 16, GMT_GOT_RADIUS = 32, }; /*! . */ enum gmt_ends { BEG = 0, END = 1 }; /*! Internal struct used in the processing of CPT z-scaling and truncation */ struct CPT_Z_SCALE { unsigned int z_adjust; /* 1 if +u was parsed and scale set, 3 if z has been adjusted, 0 otherwise */ unsigned int z_mode; /* 1 if +U was parsed, 0 otherwise */ unsigned int z_unit; /* Unit enum specified via +u */ double z_unit_to_meter; /* Scale, given z_unit, to convert z from to meters */ }; static char *GMT_just_code[12] = {"--", "LB", "CB", "RB", "--", "LM", "CM", "RM", "--", "LT", "CT", "RT"}; #define gmt_M_uneven_interval(unit) ((unit == 'o' || unit == 'O' || unit == 'k' || unit == 'K' || unit == 'R' || unit == 'r' || unit == 'D' || unit == 'd') ? true : false) /* true for uneven units */ /** @brief XYZ color of the D65 white point */ #define WHITEPOINT_X 0.950456 #define WHITEPOINT_Y 1.0 #define WHITEPOINT_Z 1.088754 /** * @brief sRGB gamma correction, transforms R to R' * http://en.wikipedia.org/wiki/SRGB */ #define GAMMACORRECTION(t) (((t) <= 0.0031306684425005883) ? (12.92*(t)) : (1.055*pow((t), 0.416666666666666667) - 0.055)) /** * @brief Inverse sRGB gamma correction, transforms R' to R */ #define INVGAMMACORRECTION(t) (((t) <= 0.0404482362771076) ? ((t)/12.92) : pow(((t) + 0.055)/1.055, 2.4)) /** * @brief CIE L*a*b* f function (used to convert XYZ to L*a*b*) * http://en.wikipedia.org/wiki/Lab_color_space */ #define LABF(t) ((t >= 8.85645167903563082e-3) ? pow(t,0.333333333333333) : (841.0/108.0)*(t) + (4.0/29.0)) #define LABINVF(t) ((t >= 0.206896551724137931) ? ((t)*(t)*(t)) : (108.0/841.0)*((t) - (4.0/29.0))) static unsigned char gmt_M_color_rgb[GMT_N_COLOR_NAMES][3] = { /* r/g/b of X11 colors */ #include "gmt_color_rgb.h" }; /*! Names of pens and their thicknesses */ struct GMT_PEN_NAME { char name[16]; double width; }; /*! . */ static struct GMT_PEN_NAME GMT_penname[GMT_N_PEN_NAMES] = { /* Names and widths of pens */ #include "gmt_pennames.h" }; /*! List of GMT color tables available in makecpt and grd2cpt * Some are based on MATLAB and Matplotlib color schemes. * To add more tables, place the master CPT file in share/cpt and * add a one-line entry with name and explanation in gmt_cpt_masters.h. */ static char *GMT_CPT_master[GMT_N_CPT_MASTERS] = { #include "gmt_cpt_masters.h" }; /* Local functions needed for public functions below */ GMT_LOCAL int gmtsupport_parse_pattern_new (struct GMT_CTRL *GMT, char *line, struct GMT_FILL *fill) { /* Parse the fill pattern syntax: p|P[+r][+b|-][+f|-] */ char *c = NULL; unsigned int first = 1; double fb_rgb[4]; fill->dpi = irint (PSL_DOTS_PER_INCH_PATTERN); if ((c = strchr (line, '+'))) { /* Got modifiers */ unsigned int pos = 0, uerr = 0; char p[GMT_BUFSIZ] = {""}; while (gmt_getmodopt (GMT, 0, c, "bfr", &pos, p, &uerr) && uerr == 0) { /* Looking for +b, +f, +r */ switch (p[0]) { case 'b': /* Background color. Giving no argument means transparent [also checking for obsolete -] */ if (p[1] == '\0' || p[1] == '-') { /* Transparent */ fill->b_rgb[0] = fill->b_rgb[1] = fill->b_rgb[2] = -1, fill->b_rgb[3] = 0; GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Background pixels set to transparent!\n"); } else { if (gmt_getrgb (GMT, &p[1], fill->b_rgb)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Background colorizing value %s not recognized!\n", &p[1]); return 2; } GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Background pixels set to colors %s\n", gmt_putrgb (GMT, fill->b_rgb)); } break; case 'f': /* Foreround color. Giving no argument means transparent [also checking for obsolete -] */ if (p[1] == '\0' || p[1] == '-') { /* Transparent */ fill->f_rgb[0] = fill->f_rgb[1] = fill->f_rgb[2] = -1, fill->f_rgb[3] = 0; GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Foreground pixels set to transparent!\n"); } else { if (gmt_getrgb (GMT, &p[1], fill->f_rgb)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Foreground colorizing value %s not recognized!\n", &p[1]); return 2; } GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Foreground pixels set to colors %s\n", gmt_putrgb (GMT, fill->f_rgb)); } break; case 'r': /* Dots-per-inch resolution */ if (p[1] == '-') { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pattern dpi %s is negative!\n", &p[1]); return 4; } fill->dpi = atoi (&p[1]); GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Pattern dpi set to %d\n", fill->dpi); break; default: break; } } if (uerr) return (GMT_PARSE_ERROR); } /* Copy name (or number) of pattern */ if (c) c[0] = '\0'; /* Chop off the modifiers */ if (!gmt_M_file_is_memory (&line[1]) && line[1] == '@') { /* Must be a cache file */ first = gmt_download_file_if_not_found (GMT, &line[1], 0) + 1; /* Add one since we started at 1 */ } strncpy (fill->pattern, &line[first], PATH_MAX-1); /* Attempt to convert to integer - will be 0 if not an integer and then we set it to -1 for a filename */ fill->pattern_no = atoi (fill->pattern); if (fill->pattern_no == 0) { fill->pattern_no = -1; gmt_set_pad (GMT, 0); /* No padding */ GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Pattern image is in file %s\n", fill->pattern); if ((fill->I = GMT_Read_Data (GMT->parent, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, fill->pattern, NULL)) == NULL) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to read image %s, no pattern set\n", fill->pattern); return (GMT_RUNTIME_ERROR); } gmt_set_pad (GMT, GMT->parent->pad); /* Restore to GMT Defaults */ fill->dim[0] = fill->I->header->n_columns; fill->dim[1] = fill->I->header->n_rows; if (fill->I->colormap) { /* Got indexed color image, must build rgb stream instead */ /* Convert colormap from integer to unsigned char and count colors */ unsigned char *colormap = gmt_M_memory (GMT, NULL, 4*256, unsigned char); int64_t n, j, k; for (n = 0; n < 4 * 256 && fill->I->colormap[n] >= 0; n++) colormap[n] = (unsigned char)fill->I->colormap[n]; n /= 4; /* Expand 8-bit indexed image to 24-bit image */ fill->I->data = gmt_M_memory (GMT, fill->I->data, 3 * fill->I->header->nm, unsigned char); n = 3 * fill->I->header->nm - 1; for (j = (int64_t)fill->I->header->nm - 1; j >= 0; j--) { k = 4 * fill->I->data[j] + 3; fill->I->data[n--] = colormap[--k], fill->I->data[n--] = colormap[--k], fill->I->data[n--] = colormap[--k]; } gmt_M_free (GMT, colormap); fill->I->header->n_bands = 3; fill->I->header->size *= 3; } else if (fill->I->header->n_bands == 4) { /* RGBA image, with a color map */ uint64_t n4, j4; for (j4 = n4 = 0; j4 < 4 * fill->I->header->nm; j4++) /* Reduce image from 32- to 24-bit */ fill->I->data[n4++] = fill->I->data[j4++], fill->I->data[n4++] = fill->I->data[j4++], fill->I->data[n4++] = fill->I->data[j4++]; fill->I->header->n_bands = 3; } fill->image = fill->I->data; fill->dim[2] = (fill->I->header->n_bands == 3) ? 24 : 8; } else GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Pattern number %d selected\n", fill->pattern_no); /* If inverse, simply flip the colors around */ if (line[0] == 'p') { GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Pattern will be inverted\n"); gmt_M_rgb_copy (fb_rgb, fill->f_rgb); gmt_M_rgb_copy (fill->f_rgb, fill->b_rgb); gmt_M_rgb_copy (fill->b_rgb, fb_rgb); } if (c) c[0] = '+'; /* Undo previous damage */ return (0); } GMT_LOCAL int gmtsupport_parse_pattern_old (struct GMT_CTRL *GMT, char *line, struct GMT_FILL *fill) { /* Parse the old-style pattern syntax */ int n, i, len, pos, end; unsigned int first = 1; char f, word[GMT_LEN256] = {""}; double fb_rgb[4]; /* We only get here if there are no +r,+f,+b strings or : in line, and we know line[1] is integer. * However, if user gave -Gp1image.jpg then we will fail trying to parse as old. So we do a check * here if the argument can be accessed as a file and if so we call the new parser and return */ if (!gmt_M_file_is_memory (&line[1]) && line[1] == '@') { /* Must be a cache file */ first = gmt_download_file_if_not_found (GMT, &line[1], 0) + 1; /* Add one since we started at 1 */ } if (!gmt_access (GMT, &line[first], F_OK)) return (gmtsupport_parse_pattern_new (GMT, line, fill)); n = sscanf (&line[1], "%d/%s", &fill->dpi, fill->pattern); if (n != 2) return (1); if (!gmt_M_file_is_memory (fill->pattern) && fill->pattern[0] == '@') { /* Must be a cache file */ first = gmt_download_file_if_not_found (GMT, fill->pattern, 0); /* Shuffle pattern name to skip the leading @ */ for (i = 1; fill->pattern[i]; i++) fill->pattern[i-1] = fill->pattern[i]; fill->pattern[i-1] = '\0'; } /* Determine if there are colorizing options applied, i.e. [:FB] */ len = (int)MIN(strlen (fill->pattern),PATH_MAX) - 1; for (i = 0, pos = -1; i < len && fill->pattern[i] && pos == -1; i++) if (i < len && fill->pattern[i] == ':' && (fill->pattern[i+1] == 'B' || fill->pattern[i+1] == 'F')) pos = i; /* THe extra i < len is needed to defeat cppcheck confusion */ if (pos > -1) fill->pattern[pos] = '\0'; fill->pattern_no = atoi (fill->pattern); if (fill->pattern_no == 0) fill->pattern_no = -1; /* See if fore- and background colors are given */ len = (int)strlen (line); for (i = 0, pos = -1; line[i] && pos == -1; i++) if (line[i] == ':' && i < len && (line[i+1] == 'B' || line[i+1] == 'F')) pos = i; pos++; if (pos > 0 && line[pos]) { /* Gave colors */ while (line[pos]) { f = line[pos++]; if (line[pos] == '-') /* Signal for transparency masking */ fb_rgb[0] = fb_rgb[1] = fb_rgb[2] = -1, fb_rgb[3] = 0; else { end = pos; while (line[end] && !(line[end] == 'F' || line[end] == 'B')) end++; strncpy (word, &line[pos], (size_t)(end - pos)); word[end - pos] = '\0'; if (gmt_getrgb (GMT, word, fb_rgb)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Colorizing value %s not recognized!\n", word); return 2; } } if (f == 'f' || f == 'F') gmt_M_rgb_copy (fill->f_rgb, fb_rgb); else if (f == 'b' || f == 'B') gmt_M_rgb_copy (fill->b_rgb, fb_rgb); else { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Colorizing argument %c not recognized!\n", f); return 2; } while (line[pos] && !(line[pos] == 'F' || line[pos] == 'B')) pos++; } } /* If inverse, simply flip the colors around */ if (line[0] == 'p') { gmt_M_rgb_copy (fb_rgb, fill->f_rgb); gmt_M_rgb_copy (fill->f_rgb, fill->b_rgb); gmt_M_rgb_copy (fill->b_rgb, fb_rgb); } return 0; } /*! . */ GMT_LOCAL int gmtsupport_parse_pattern (struct GMT_CTRL *GMT, char *line, struct GMT_FILL *fill) { int err; /* New syntax may have a modifier */ if (strstr(line, "+r") || strstr(line, "+f") || strstr(line, "+b") || !strchr(line, '/')) /* Clearly new syntax */ err = gmtsupport_parse_pattern_new (GMT, line, fill); else err = gmtsupport_parse_pattern_old (GMT, line, fill); fill->use_pattern = true; return (err); } /*! . */ GMT_LOCAL char *gmtsupport_get_userimagename (struct GMT_CTRL *GMT, char *line, char *cpt_path) { /* When a cpt is not in the current directory but given by relative or absolute path * AND that cpt refers to a user pattern file (which may be a relative or absolute path) * then unless that pattern file can be found we will try to prepend the path to the cpt * file and see if the pattern can be found that way. */ int j, err; char *name = NULL, path[PATH_MAX+GMT_LEN256] = {""}; struct GMT_FILL fill; if (!gmt_M_is_pattern (line)) return NULL; /* Not an image specification */ err = gmtsupport_parse_pattern (GMT, line, &fill); /* See if this returns an error or not */ if (err) return NULL; /* Not a valid image specification */ if (fill.pattern_no > 0) return NULL; /* Not a user image */ /* Here we do have a pattern specification */ /* Try the user's default directories */ if (gmtlib_getuserpath (GMT, fill.pattern, path)) return NULL; /* Yes, found so no problems */ /* Now must put our faith in the cpt path and hope it has a path that can help us */ if (cpt_path == NULL || cpt_path[0] == '<') { /* Without an actual file path we must warn and bail */ GMT_Report (GMT->parent, GMT_MSG_ERROR, "Not enough information to determine location of user pattern %s\n", fill.pattern); return NULL; } j = (int)strlen (cpt_path); while (j > 0 && cpt_path[j] != '/') j--; /* Find last slash */ if (j > 0) { /* OK, got the cpt directory */ cpt_path[j] = '\0'; /* Temporarily chop off the slash */ sprintf (path, "%s/%s", cpt_path, fill.pattern); cpt_path[j] = '/'; /* Restore the slash */ if (access (path, R_OK)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Not enough information to determine location of user pattern %s\n", fill.pattern); return NULL; } name = strdup (path); /* Must use this image name instead as it contains the full working path */ } return (name); } /*! . */ void gmt_rgb_to_hsv (double rgb[], double hsv[]) { double diff; unsigned int i, imax = 0, imin = 0; /* This had checks using rgb value in doubles (e.g. (max_v == xr)), which failed always on some compilers. Changed to integer logic: 2009-02-05 by RS. */ hsv[3] = rgb[3]; /* Pass transparency unchanged */ for (i = 1; i < 3; i++) { if (rgb[i] > rgb[imax]) imax = i; if (rgb[i] < rgb[imin]) imin = i; } diff = rgb[imax] - rgb[imin]; hsv[0] = 0.0; hsv[1] = (rgb[imax] == 0.0) ? 0.0 : diff / rgb[imax]; hsv[2] = rgb[imax]; if (hsv[1] == 0.0) return; /* Hue is undefined */ hsv[0] = 120.0 * imax + 60.0 * (rgb[(imax + 1) % 3] - rgb[(imax + 2) % 3]) / diff; if (hsv[0] < 0.0) hsv[0] += 360.0; if (hsv[0] > 360.0) hsv[0] -= 360.0; } /*! . */ void gmt_hsv_to_rgb (double rgb[], double hsv[]) { int i; double h, f, p, q, t, rr, gg, bb; rgb[3] = hsv[3]; /* Pass transparency unchanged */ if (hsv[1] == 0.0) rgb[0] = rgb[1] = rgb[2] = hsv[2]; else { h = hsv[0]; while (h >= 360.0) h -= 360.0; while (h < 0.0) h += 360.0; h /= 60.0; i = irint (floor (h)); f = h - i; p = hsv[2] * (1.0 - hsv[1]); q = hsv[2] * (1.0 - (hsv[1] * f)); t = hsv[2] * (1.0 - (hsv[1] * (1.0 - f))); switch (i) { case 0: rr = hsv[2]; gg = t; bb = p; break; case 1: rr = q; gg = hsv[2]; bb = p; break; case 2: rr = p; gg = hsv[2]; bb = t; break; case 3: rr = p; gg = q; bb = hsv[2]; break; case 4: rr = t; gg = p; bb = hsv[2]; break; default: rr = hsv[2]; gg = p; bb = q; break; } rgb[0] = (rr < 0.0) ? 0.0 : rr; rgb[1] = (gg < 0.0) ? 0.0 : gg; rgb[2] = (bb < 0.0) ? 0.0 : bb; } } /*! . */ GMT_LOCAL void gmtsupport_rgb_to_cmyk (double rgb[], double cmyk[]) { /* Plain conversion; with default undercolor removal or blackgeneration */ /* RGB is in 0-1, CMYK will be in 0-1 range */ int i; cmyk[4] = rgb[3]; /* Pass transparency unchanged */ for (i = 0; i < 3; i++) cmyk[i] = 1.0 - rgb[i]; cmyk[3] = MIN (cmyk[0], MIN (cmyk[1], cmyk[2])); /* Default Black generation */ if (cmyk[3] < GMT_CONV8_LIMIT) cmyk[3] = 0.0; /* To implement device-specific blackgeneration, supply lookup table K = BG[cmyk[3]] */ for (i = 0; i < 3; i++) { cmyk[i] -= cmyk[3]; /* Default undercolor removal */ if (cmyk[i] < GMT_CONV8_LIMIT) cmyk[i] = 0.0; } /* To implement device-specific undercolor removal, supply lookup table u = UR[cmyk[3]] */ } /*! . */ GMT_LOCAL void gmtsupport_cmyk_to_rgb (double rgb[], double cmyk[]) { /* Plain conversion; no undercolor removal or blackgeneration */ /* CMYK is in 0-1, RGB will be in 0-1 range */ int i; rgb[3] = cmyk[4]; /* Pass transparency unchanged */ for (i = 0; i < 3; i++) rgb[i] = 1.0 - cmyk[i] - cmyk[3]; } /*! . */ GMT_LOCAL void gmtsupport_cmyk_to_hsv (double hsv[], double cmyk[]) { /* Plain conversion; no undercolor removal or blackgeneration */ /* CMYK is in 0-1, RGB will be in 0-1 range */ double rgb[4]; gmtsupport_cmyk_to_rgb (rgb, cmyk); gmt_rgb_to_hsv (rgb, hsv); } /*! * Transform sRGB to CIE XYZ with the D65 white point * * Poynton, "Frequently Asked Questions About Color," page 10 * Wikipedia: http://en.wikipedia.org/wiki/SRGB * Wikipedia: http://en.wikipedia.org/wiki/CIE_1931_color_space */ void gmt_rgb_to_xyz (double rgb[], double xyz[]) { double R, G, B; R = INVGAMMACORRECTION(rgb[0]); G = INVGAMMACORRECTION(rgb[1]); B = INVGAMMACORRECTION(rgb[2]); xyz[0] = 0.4123955889674142161*R + 0.3575834307637148171*G + 0.1804926473817015735*B; xyz[1] = 0.2125862307855955516*R + 0.7151703037034108499*G + 0.07220049864333622685*B; xyz[2] = 0.01929721549174694484*R + 0.1191838645808485318*G + 0.9504971251315797660*B; } /*! . */ void gmt_xyz_to_rgb (double rgb[], double xyz[]) { double R, G, B, min; R = ( 3.2406 * xyz[0] - 1.5372 * xyz[1] - 0.4986 * xyz[2]); G = (-0.9689 * xyz[0] + 1.8758 * xyz[1] + 0.0415 * xyz[2]); B = ( 0.0557 * xyz[0] - 0.2040 * xyz[1] + 1.0570 * xyz[2]); min = MIN(MIN(R, G), B); /* Force nonnegative values so that gamma correction is well-defined. */ if (min < 0) { R -= min; G -= min; B -= min; } /* Transform from RGB to R'G'B' */ rgb[0] = GAMMACORRECTION(R); rgb[1] = GAMMACORRECTION(G); rgb[2] = GAMMACORRECTION(B); } /*! * Convert CIE XYZ to CIE L*a*b* (CIELAB) with the D65 white point * * Wikipedia: http://en.wikipedia.org/wiki/Lab_color_space */ void gmt_xyz_to_lab (double xyz[], double lab[]) { double X, Y, Z; X = LABF( xyz[0] / WHITEPOINT_X ); Y = LABF( xyz[1] / WHITEPOINT_Y ); Z = LABF( xyz[2] / WHITEPOINT_Z ); lab[0] = 116 * Y - 16; lab[1] = 500 * (X - Y); lab[2] = 200 * (Y - Z); } /*! . */ void gmt_lab_to_xyz (double xyz[], double lab[]) { xyz[0] = WHITEPOINT_X * LABINVF( (lab[0] + 16)/116 + lab[1]/500 ); xyz[1] = WHITEPOINT_Y * LABINVF( (lab[0] + 16)/116 ); xyz[2] = WHITEPOINT_Z * LABINVF( (lab[0] + 16)/116 - lab[2]/200 ); } /*! . */ void gmt_rgb_to_lab (double rgb[], double lab[]) { /* RGB is in 0-1, LAB will be in ??? range */ double xyz[3]; gmt_rgb_to_xyz (rgb, xyz); gmt_xyz_to_lab (xyz, lab); } /*! . */ void gmt_lab_to_rgb (double rgb[], double lab[]) { double xyz[3]; gmt_lab_to_xyz (xyz, lab); gmt_xyz_to_rgb (rgb, xyz); } /*! . */ GMT_LOCAL int gmtsupport_comp_double_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ bool bad_1, bad_2; const double *point_1 = p_1, *point_2 = p_2; bad_1 = gmt_M_is_dnan ((*point_1)); bad_2 = gmt_M_is_dnan ((*point_2)); if (bad_1 && bad_2) return (0); if (bad_1) return (+1); if (bad_2) return (-1); if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_float_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ bool bad_1, bad_2; const float *point_1 = p_1, *point_2 = p_2; bad_1 = gmt_M_is_fnan ((*point_1)); bad_2 = gmt_M_is_fnan ((*point_2)); if (bad_1 && bad_2) return (0); if (bad_1) return (+1); if (bad_2) return (-1); if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_ulong_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const uint64_t *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_long_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const int64_t *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_uint_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const unsigned int *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_int_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const int *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_ushort_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const unsigned short int *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_short_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const short int *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_uchar_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const unsigned char *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL int gmtsupport_comp_char_asc (const void *p_1, const void *p_2) { /* Returns -1 if point_1 is < that point_2, +1 if point_2 > point_1, and 0 if they are equal */ const char *point_1 = p_1, *point_2 = p_2; if ((*point_1) < (*point_2)) return (-1); if ((*point_1) > (*point_2)) return (+1); return (0); } /*! . */ GMT_LOCAL bool gmtsupport_check_irgb (int irgb[], double rgb[]) { if ((irgb[0] < 0 || irgb[0] > 255) || (irgb[1] < 0 || irgb[1] > 255) || (irgb[2] < 0 || irgb[2] > 255)) return (true); rgb[0] = gmt_M_is255 (irgb[0]); rgb[1] = gmt_M_is255 (irgb[1]); rgb[2] = gmt_M_is255 (irgb[2]); return (false); } /*! . */ GMT_LOCAL bool gmtsupport_check_rgb (double rgb[]) { return ((rgb[0] < 0.0 || rgb[0] > 1.0) || (rgb[1] < 0.0 || rgb[1] > 1.0) || (rgb[2] < 0.0 || rgb[2] > 1.0)); } /*! . */ GMT_LOCAL bool gmtsupport_check_hsv (double hsv[]) { return ((hsv[0] < 0.0 || hsv[0] > 360.0) || (hsv[1] < 0.0 || hsv[1] > 1.0) || (hsv[2] < 0.0 || hsv[2] > 1.0)); } /*! . */ GMT_LOCAL bool gmtsupport_check_cmyk (double cmyk[]) { unsigned int i; for (i = 0; i < 4; i++) cmyk[i] *= 0.01; for (i = 0; i < 4; i++) if (cmyk[i] < 0.0 || cmyk[i] > 1.0) return (true); return (false); } /*! . */ unsigned int gmt_char_count (char *txt, char c) { unsigned int i = 0, n = 0; while (txt[i]) if (txt[i++] == c) n++; return (n); } /*! . */ GMT_LOCAL bool gmtsupport_gethsv (struct GMT_CTRL *GMT, char *line, double hsv[]) { int n, i, count, irgb[3], c = 0; double rgb[4], cmyk[5]; char buffer[GMT_LEN64] = {""}, *t = NULL; if (!line) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "No argument given to gmtsupport_gethsv\n"); GMT->parent->error = GMT_PARSE_ERROR; return false; } if (!line[0]) return (false); /* Nothing to do - accept default action */ rgb[3] = hsv[3] = cmyk[4] = 0.0; /* Default is no transparency */ if (line[0] == '-') { hsv[0] = -1.0; hsv[1] = -1.0; hsv[2] = -1.0; return (false); } strncpy (buffer, line, GMT_LEN64-1); /* Make local copy */ if ((t = strstr (buffer, "@")) && strlen (t) > 1) { /* User requested transparency via @ */ double transparency = atof (&t[1]); if (transparency < 0.0 || transparency > 100.0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "Representation of transparency (%s) not recognized. Using default [0 or opaque].\n", t); else rgb[3] = hsv[3] = cmyk[4] = transparency / 100.0; /* Transparency is in 0-1 range */ t[0] = '\0'; /* Chop off transparency for the rest of this function */ } if (buffer[0] == '#') { /* #rrggbb */ n = sscanf (buffer, "#%2x%2x%2x", (unsigned int *)&irgb[0], (unsigned int *)&irgb[1], (unsigned int *)&irgb[2]); if (n != 3 || gmtsupport_check_irgb (irgb, rgb)) return (true); gmt_rgb_to_hsv (rgb, hsv); return (false); } /* If it starts with a letter, then it could be a name */ if (isalpha ((unsigned char) buffer[0])) { if ((n = (int)gmt_colorname2index (GMT, buffer)) < 0) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Colorname %s not recognized!\n", buffer); return (true); } for (i = 0; i < 3; i++) rgb[i] = gmt_M_is255 (gmt_M_color_rgb[n][i]); gmt_rgb_to_hsv (rgb, hsv); return (false); } /* Definitely wrong, at this point, is something that does not end in a number */ if (strlen(buffer) < 1) return (true); /* Nothing, which is bad */ c = buffer[strlen(buffer)-1]; if (!(isdigit (c) || c == '.')) return (true); count = (int)gmt_char_count (buffer, '/'); if (count == 3) { /* c/m/y/k */ n = sscanf (buffer, "%lf/%lf/%lf/%lf", &cmyk[0], &cmyk[1], &cmyk[2], &cmyk[3]); if (n != 4 || gmtsupport_check_cmyk (cmyk)) return (true); gmtsupport_cmyk_to_hsv (hsv, cmyk); return (false); } if (count == 2) { /* r/g/b */ n = sscanf (buffer, "%d/%d/%d", &irgb[0], &irgb[1], &irgb[2]); if (n != 3 || gmtsupport_check_irgb (irgb, rgb)) return (true); gmt_rgb_to_hsv (rgb, hsv); return (false); } if (gmt_char_count (buffer, '-') == 2) { /* h-s-v */ n = sscanf (buffer, "%lf-%lf-%lf", &hsv[0], &hsv[1], &hsv[2]); return (n != 3 || gmtsupport_check_hsv (hsv)); } if (count == 0) { /* gray */ n = sscanf (buffer, "%d", &irgb[0]); irgb[1] = irgb[2] = irgb[0]; if (n != 1 || gmtsupport_check_irgb (irgb, rgb)) return (true); gmt_rgb_to_hsv (rgb, hsv); return (false); } /* Get here if there is a problem */ return (true); } /*! . */ GMT_LOCAL bool gmtsupport_is_pattern_old (char *word) { /* Returns true if the word is a pattern specification * Old style: P|p/[:B[F]] */ if (strchr (word, ':')) return (true); /* Only patterns may have a colon */ if (!(word[0] == 'P' || word[0] == 'p')) return false; /* Patterns must start with P or p */ if (!strchr (word, '/')) return (false); /* Patterns separate dpi and pattern with a slash */ /* Here we know we start with P|p and there is a slash - this can only be a pattern specification */ return (true); } /*! . */ GMT_LOCAL bool gmtsupport_is_pattern_new (struct GMT_CTRL *GMT, char *word) { /* Returns true if the word is a pattern specification * New style: P|p[+b][+f][+r] */ char *c = NULL; int n; if ((c = strchr (word, '+')) && strchr ("bfr", c[1])) return (true); /* Found +b, +f or +r */ /* Here we have no modifiers at all */ if (!(word[0] == 'P' || word[0] == 'p')) return false; /* Patterns must start with P or p */ /* Pattern without modifier must be an integer OR a valid file */ if ((n = atoi (&word[1])) > 0 && n < PSL_N_PATTERNS) return (true); /* Got a valid integer */ /* Se if we can access the file */ if (!gmt_access (GMT, &word[1], F_OK)) return (true); /* Got a file that exists */ return (false); /* Not a pattern */ } /*! . */ GMT_LOCAL bool gmtsupport_is_pattern (struct GMT_CTRL *GMT, char *word) { /* Returns true if the word is a pattern specification * Old style: P|p/[:B[F]] * New style: P|p[+b][+f][+r] */ bool val; /* New syntax may have a modifier or no slash AND no colon */ if ((strstr(word, "+r") || strstr(word, "+f") || strstr(word, "+b") || !strchr(word, '/')) && !strchr (word,':')) val = gmtsupport_is_pattern_new (GMT, word); else val = gmtsupport_is_pattern_old (word); return (val); } /*! . */ bool gmtlib_is_color (struct GMT_CTRL *GMT, char *word) { int i, k, n, n_hyphen = 0, n_slashes = 0; /* Returns true if we are sure the word is a color string - else false. * color syntax is |||. * NOTES: 1) is excluded since this function is called in places where * a single integer may be used for font size or pen width... * 2) We are not checking if the values are kosher; just that they follow the pattern */ n = (int)strlen (word); if (n == 0) return (false); if (word[0] == '#') return (true); /* Probably #rrggbb */ if (gmt_colorname2index (GMT, word) >= 0) return (true); /* Valid color name */ /* Skip dimension specifications with units c|i|m|p or a colon */ if (strchr(word,'t')) return (false); /* Got a t somewhere */ if (strchr(word,':')) return (false); /* Got a : somewhere */ if (strchr(word,'c')) return (false); /* Got a c somewhere */ if (strchr(word,'i')) return (false); /* Got an i somewhere */ if (strchr(word,'m')) return (false); /* Got a m somewhere */ if (strchr(word,'p')) return (false); /* Got a p somewhere */ for (i = k = 0; word[i]; i++) if (word[i] == '/') n_slashes++; /* Count slashes */ if (n_slashes == 1 || n_slashes > 3) return (false); /* No color spec takes only 1 slash or more than 3 */ n--; /* Check for h-s-v as well. Must find 2 hyphens */ while (n >= 0 && (strchr ("/-.", word[n]) || isdigit ((int)word[n]))) { if (word[n] == '-') n_hyphen++; n--; /* Wind down as long as we find -,.,/ or digits */ } return ((n == -1 && (n_slashes || n_hyphen == 2))); /* Hence will fail the slash test */ } /*! . */ int gmt_getfonttype (struct GMT_CTRL *GMT, char *name) { unsigned int i; if (!name[0]) return (-1); if (!isdigit ((unsigned char) name[0])) { /* Does not start with number. Try font name */ int ret = -1; for (i = 0; i < GMT->session.n_fonts && strcmp (name, GMT->session.font[i].name); i++); if (i < GMT->session.n_fonts) ret = i; return (ret); } if (!isdigit ((unsigned char) name[strlen(name)-1])) return (-1); /* Starts with digit, ends with something else: cannot be */ return (atoi (name)); } /*! . */ GMT_LOCAL bool gmtsupport_is_fontname (struct GMT_CTRL *GMT, char *word) { /* Returns true if the word is one of the named fonts */ unsigned int i; if (!word[0]) return (false); for (i = 0; i < GMT->session.n_fonts && strcmp (word, GMT->session.font[i].name); i++); if (i == GMT->session.n_fonts) return (false); return (true); } /*! . */ GMT_LOCAL int gmtsupport_name2pen (char *name) { /* Return index into structure with pennames and width, for given name */ int i, k; char Lname[GMT_LEN64] = {""}; strncpy (Lname, name, GMT_LEN64-1); gmt_str_tolower (Lname); for (i = 0, k = -1; k < 0 && i < GMT_N_PEN_NAMES; i++) if (!strcmp (Lname, GMT_penname[i].name)) k = i; return (k); } /*! . */ GMT_LOCAL int gmtsupport_pen2name (double width) { /* Return index into structure with pennames and width, for given width */ int i, k; for (i = 0, k = -1; k < 0 && i < GMT_N_PEN_NAMES; i++) if (gmt_M_eq (width, GMT_penname[i].width)) k = i; return (k); } /*! . */ GMT_LOCAL int gmtsupport_getpenwidth (struct GMT_CTRL *GMT, char *line, struct GMT_PEN *P) { int n; /* SYNTAX for pen width: [p|i|c|m] or [fat, thin, etc] */ if (!line || !line[0]) return (GMT_NOERROR); /* Nothing given, return */ if ((line[0] == '.' && (line[1] >= '0' && line[1] <= '9')) || (line[0] >= '0' && line[0] <= '9')) { /* Pen thickness with optional unit at end */ P->width = gmt_convert_units (GMT, line, GMT_PT, GMT_PT); } else { /* Pen name was given - these refer to fixed widths in points */ if ((n = gmtsupport_name2pen (line)) < 0) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pen name %s not recognized!\n", line); return GMT_PARSE_ERROR; } P->width = GMT_penname[n].width; } return (GMT_NOERROR); } /*! . */ GMT_LOCAL int gmtsupport_getpenstyle (struct GMT_CTRL *GMT, char *line, struct GMT_PEN *P) { unsigned int i, n, pos, unit = GMT_PT, n_dash = 0; double width; char tmp[GMT_PEN_LEN] = {""}, string[GMT_BUFSIZ] = {""}, ptr[GMT_BUFSIZ] = {""}; if (!line || !line[0]) return (GMT_NOERROR); /* Nothing to do as no style given. Note we do not reset any existing style here; use solid to do so */ if (!strcmp (line, "solid")) { /* Used to override any other default style */ P->offset = 0.0; P->style[0] = '\0'; return (GMT_NOERROR); } if (!strncmp (line, "dashdot", 7U)) strcpy (line, "-."); /* Accept "dashdot*" to mean -. */ if (!strncmp (line, "dotdash", 7U)) strcpy (line, ".-"); /* Accept "dotdash*" to mean .- */ if (!strncmp (line, "dash", 4U)) strcpy (line, "-"); /* Accept "dash*" to mean - */ if (!strncmp (line, "dot", 3U)) strcpy (line, "."); /* Accept "dot*" to mean . */ if (!strcmp (line, "a")) { /* Old GMT4 "a" style */ if (gmt_M_compat_check (GMT, 4)) GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Pen-style \"a\" is deprecated, use \"dashed\" or \"-\" instead\n"); else { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad pen-style %s\n", line); return GMT_PARSE_ERROR; } strcpy (line, "-"); /* Accepted GMT4 style "a" to mean - */ } else if (!strcmp (line, "o")) { /* Old GMT4 "o" style */ if (gmt_M_compat_check (GMT, 4)) GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Pen-style \"o\" is deprecated, use \"dotted\" or \".\" instead\n"); else { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad pen-style %s\n", line); return GMT_PARSE_ERROR; } strcpy (line, "."); /* Accepted GMT4 style "a" to mean - */ } n = (int)strlen (line) - 1; if (strchr (GMT_DIM_UNITS, line[n])) /* Separate unit given to style string */ unit = gmtlib_unit_lookup (GMT, line[n], GMT->current.setting.proj_length_unit); width = (P->width < GMT_CONV4_LIMIT) ? GMT_PENWIDTH : P->width; if (isdigit ((int)line[0])) { /* Specified numeric pattern will start with an integer */ unsigned int c_pos; for (i = 1, c_pos = 0; line[i] && c_pos == 0; i++) if (line[i] == ':') c_pos = i; if (c_pos == 0) { GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Pen style %s does not contain :. set to 0\n", line); P->offset = 0.0; sscanf (line, "%s", P->style); } else { line[c_pos] = ' '; sscanf (line, "%s %lf", P->style, &P->offset); line[c_pos] = ':'; } for (i = 0; P->style[i]; i++) if (P->style[i] == '_') P->style[i] = ' '; /* Must convert given units to points */ pos = 0; while ((gmt_strtok (P->style, " ", &pos, ptr))) { snprintf (tmp, GMT_PEN_LEN, "%g ", (atof (ptr) * GMT->session.u2u[unit][GMT_PT])); strcat (string, tmp); n_dash++; } string[strlen (string) - 1] = 0; if (strlen (string) >= GMT_PEN_LEN) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pen attributes too long!\n"); return GMT_PARSE_ERROR; } strncpy (P->style, string, GMT_PEN_LEN-1); P->offset *= GMT->session.u2u[unit][GMT_PT]; } else { /* New way of building it up with - and . */ P->style[0] = '\0'; P->offset = 0.0; for (i = 0; line[i]; i++) { if (line[i] == '-') { /* Dash */ snprintf (tmp, GMT_PEN_LEN, "%g %g ", 8.0 * width, 4.0 * width); strcat (P->style, tmp); n_dash += 2; } else if (line[i] == '.') { /* Dot */ snprintf (tmp, GMT_PEN_LEN, "%g %g ", width, 4.0 * width); strcat (P->style, tmp); n_dash += 2; } else { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pen attributes not using just - and . for dashes and dots. Offending character --> %c\n", line[i]); return GMT_PARSE_ERROR; } } P->style[strlen(P->style)-1] = '\0'; /* Chop off trailing space */ } if (n_dash >= 11) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pen attributes contain more than 11 items (limit for PostScript setdash operator)!\n"); return GMT_PARSE_ERROR; } return (GMT_NOERROR); } /*! . */ GMT_LOCAL bool gmtsupport_is_penstyle (char *word) { int n; /* Returns true if we are sure the word is a style string - else false. * style syntax is a|o|[:]| * Also recognized "dashed" for -, "dotted" for . as well as "solid" */ n = (int)strlen (word); if (n == 0) return (false); if (!strncmp (word, "dotdash", 7U) || !strncmp (word, "dashdot", 7U) \ || !strncmp (word, "dash", 4U) || !strncmp (word, "dot", 3U) \ || !strncmp (word, "solid", 5U)) return (true); n--; if (strchr (GMT_DIM_UNITS, word[n])) n--; /* Reduce length by 1; the unit character */ if (n < 0) return (false); /* word only contained a unit character? */ if (n == 0) { if (word[0] == '-' || word[0] == 'a' || word[0] == '.' || word[0] == 'o') return (true); return (false); /* No other 1-char style patterns possible */ } if (strchr(word,'t')) return (false); /* Got a t somewhere */ if (strchr(word,':')) return (true); /* Got : */ if (strchr(word,'_')) { /* Got without optional :, add :0 */ strcat (word, ":0"); return (true); } while (n >= 0 && (word[n] == '-' || word[n] == '.')) n--; /* Wind down as long as we find - or . */ return ((n == -1)); /* true if we only found -/., false otherwise */ } /*! . */ GMT_LOCAL void gmtsupport_free_range (struct GMT_CTRL *GMT, struct GMT_LUT *S) { gmt_M_free (GMT, S->label); gmt_M_free (GMT, S->fill); } /*! . */ GMT_LOCAL void gmtsupport_copy_palette_hdrs (struct GMT_CTRL *GMT, struct GMT_PALETTE *P_to, struct GMT_PALETTE *P_from) { unsigned int hdr; if (P_from->n_headers == 0) return; /* Nothing to do */ /* Must duplicate the header records */ P_to->n_headers = P_from->n_headers; if (P_to->n_headers) P_to->header = gmt_M_memory (GMT, NULL, P_from->n_headers, char *); for (hdr = 0; hdr < P_from->n_headers; hdr++) P_to->header[hdr] = strdup (P_from->header[hdr]); } #if 0 /*! Decode the optional +w, +u|U and determine scales */ GMT_LOCAL struct CPT_Z_SCALE *gmtsupport_cpt_parse (struct GMT_CTRL *GMT, char *file, unsigned int direction) { unsigned int pos = 0, uerr = 0; int unit; char *c = NULL, p[GMT_LEN32] = {""}; struct CPT_Z_SCALE *Z = NULL; gmt_M_unused(direction); if ((c = gmt_first_modifier (GMT, file, "uUw")) == NULL) return NULL; /* No modifiers found */ /* Found at least one modifier */ Z = gmt_M_memory (GMT, NULL, 1, struct CPT_Z_SCALE); while (gmt_getmodopt (GMT, 0, c, "uUw", &pos, p, &uerr) && uerr == 0) { switch (p[0]) { case 'U': Z->z_mode = 1; /* Intentionally fall through */ case 'u': unit = gmtlib_get_unit_number (GMT, p[1]); /* Convert char unit to enumeration constant for this unit */ if (unit == GMT_IS_NOUNIT) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "CPT z unit specification %s was unrecognized (part of file name?) and is ignored.\n", c); gmt_M_free (GMT, Z); return NULL; } Z->z_unit = (unsigned int)unit; Z->z_unit_to_meter = GMT->current.proj.m_per_unit[Z->z_unit]; /* Converts unit to meters */ if (Z->z_mode) Z->z_unit_to_meter = 1.0 / Z->z_unit_to_meter; /* Wanted the inverse */ Z->z_adjust |= 1; /* Says we have successfully parsed and readied the x/y scaling */ break; case 'w': if (direction == GMT_IN) Z->z_wrap = 1; else GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmtsupport_cpt_parse: Modifier +w ignored for output.\n"); break; default: /* These are caught in gmt_getmodopt so break is just for Coverity */ break; } } c[0] = '\0'; /* Chop off all modifiers so the filename can be used */ return (Z); } #endif /*! Decode the optional +u|U and determine scales and +h[] for hinge control */ GMT_LOCAL struct CPT_Z_SCALE *gmtsupport_cpt_parse (struct GMT_CTRL *GMT, char *file, unsigned int direction, unsigned int *hinge_mode, double *z_hinge) { /* CPT file arg is [+h[]][+u|U] * The +h modifier is used to turn a soft hinge in a CPT to a hard hinge at the user-selected z-value [0] * or to adjust the location of the hinge for an hard hinge CPT. */ enum gmt_enum_units u_number; unsigned int mode = 0; char *c = NULL; struct CPT_Z_SCALE *Z = NULL; gmt_M_unused(direction); *hinge_mode = 0; /* Default is no hinge modifier */ if ((c = strstr (file, "+h"))) { /* Gave hinge modifier, examine and set parameters */ if (c[2]) { /* Gave hinge value for soft hinge */ if (gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Z), gmt_scanf (GMT, &c[2], gmt_M_type (GMT, GMT_IN, GMT_Z), z_hinge), &c[2])) { GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtsupport_cpt_parse: CPT hinge modifier %s was not successfully parsed and is ignored.\n", c); } else { /* Parsed successfully, this turns a soft CPT hinge to a hard user hinge */ c[0] = '\0'; /* Chop off the hinge specification from the file name */ GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "gmtsupport_cpt_parse: CPT hard hinge was added at z = %s for file %s\n", &c[2], file); *hinge_mode = 1; } } else { /* Accept zero as hard hinge value */ *hinge_mode = 1; c[0] = '\0'; /* Chop off the hinge specification from the file name */ GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "gmtsupport_cpt_parse: CPT CPT hard hinge was added at z = 0 for file %s\n", file); } } if ((c = gmtlib_file_unitscale (file)) == NULL) return NULL; /* Did not find any +u|U modifiers */ mode = (c[1] == 'u') ? 0 : 1; u_number = gmtlib_get_unit_number (GMT, c[2]); /* Convert char unit to enumeration constant for this unit */ if (u_number == GMT_IS_NOUNIT) { GMT_Report (GMT->parent, GMT_MSG_WARNING, "CPT z unit specification %s was unrecognized (part of file name?) and is ignored.\n", c); return NULL; } /* Got a valid unit */ Z = gmt_M_memory (GMT, NULL, 1, struct CPT_Z_SCALE); Z->z_unit_to_meter = GMT->current.proj.m_per_unit[u_number]; /* Converts unit to meters */ if (mode) Z->z_unit_to_meter = 1.0 / Z->z_unit_to_meter; /* Wanted the inverse */ Z->z_unit = u_number; /* Unit ID */ Z->z_adjust |= 1; /* Says we have successfully parsed and readied the x/y scaling */ Z->z_mode = mode; c[0] = '\0'; /* Chop off the unit specification from the file name */ return (Z); } /*! . */ GMT_LOCAL int gmtsupport_find_cpt_hinge (struct GMT_CTRL *GMT, struct GMT_PALETTE *P) { /* Return the slice number where z_low' = 0 for a CPT with hinge and normalized z' = -1 to +1 */ unsigned int k; if (!P->has_hinge) return GMT_NOTSET; /* Does not have any hinge */ for (k = 0; k < P->n_colors; k++) if (doubleAlmostEqualZero (0.0, P->data[k].z_low)) { GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Found CPT hinge at z' = 0 for slice k = %u!\n", k); return (int)k; } return GMT_NOTSET; } /*! . */ GMT_LOCAL void gmtsupport_cpt_z_scale (struct GMT_CTRL *GMT, struct GMT_PALETTE *P, struct CPT_Z_SCALE *Z, unsigned int direction) { unsigned int k; double scale = 1.0; struct GMT_PALETTE_HIDDEN *PH = gmt_get_C_hidden (P); /* Apply the scaling of z as given by the header's z_* settings. * After reading a cpt it will have z in meters. * Before writing a cpt, it may have units changed back to original units * or scaled to another set of units */ if (direction == GMT_IN) { if (Z) { PH->z_adjust[GMT_IN] = Z->z_adjust; PH->z_unit[GMT_IN] = Z->z_unit; PH->z_mode[GMT_IN] = Z->z_mode; PH->z_unit_to_meter[GMT_IN] = Z->z_unit_to_meter; } if (PH->z_adjust[GMT_IN] == 0) return; /* Nothing to do */ if (PH->z_adjust[GMT_IN] & 2) return; /* Already scaled them */ scale = PH->z_unit_to_meter[GMT_IN]; /* To multiply all z-related entries in the CPT */ PH->z_adjust[GMT_IN] = 2; /* Now the cpt is ready for use and in meters */ if (PH->z_mode[GMT_IN]) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Input CPT z unit was converted from meters to %s after reading.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_IN]]); else GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Input CPT z unit was converted from %s to meters after reading.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_IN]]); } else if (direction == GMT_OUT) { /* grid x/y are assumed to be in meters */ if (Z) {PH->z_adjust[GMT_OUT] = Z->z_adjust; PH->z_unit[GMT_OUT] = Z->z_unit; PH->z_mode[GMT_OUT] = Z->z_mode; PH->z_unit_to_meter[GMT_OUT] = Z->z_unit_to_meter; } if (PH->z_adjust[GMT_OUT] & 1) { /* Was given a new unit for output */ scale = 1.0 / PH->z_unit_to_meter[GMT_OUT]; PH->z_adjust[GMT_OUT] = 2; /* Now we are ready for writing */ if (PH->z_mode[GMT_OUT]) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output CPT z unit was converted from %s to meters before writing.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_OUT]]); else GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output CPT z unit was converted from meters to %s before writing.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_OUT]]); } else if (PH->z_adjust[GMT_IN] & 2) { /* Just undo old scaling */ scale = 1.0 / PH->z_unit_to_meter[GMT_IN]; PH->z_adjust[GMT_IN] -= 2; /* Now it is back to where we started */ if (PH->z_mode[GMT_OUT]) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output CPT z unit was reverted back to %s from meters before writing.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_IN]]); else GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Output CPT z unit was reverted back from meters to %s before writing.\n", GMT->current.proj.unit_name[PH->z_unit[GMT_IN]]); } } /* If we got here we must scale the CPT's z-values */ for (k = 0; k < P->n_colors; k++) { P->data[k].z_high *= scale; P->data[k].z_low *= scale; P->data[k].i_dz /= scale; } } /*! . */ GMT_LOCAL void gmtsupport_truncate_cpt_slice (struct GMT_LUT *S, bool do_hsv, double z_cut, int side) { /* Interpolate this slice to the new low or high value, depending on side */ double f = (z_cut - S->z_low) * S->i_dz, hsv[4], rgb[4]; unsigned int k; if (do_hsv) { /* Interpolation in HSV space */ for (k = 0; k < 4; k++) hsv[k] = S->hsv_low[k] + S->hsv_diff[k] * f; gmt_hsv_to_rgb (rgb, hsv); } else { /* Interpolation in RGB space */ for (k = 0; k < 4; k++) rgb[k] = S->rgb_low[k] + S->rgb_diff[k] * f; gmt_rgb_to_hsv (rgb, hsv); } if (side == -1) { gmt_M_memcpy (S->hsv_low, hsv, 4, double); gmt_M_memcpy (S->rgb_low, rgb, 4, double); S->z_low = z_cut; } else { /* Last slice */ gmt_M_memcpy (S->hsv_high, hsv, 4, double); gmt_M_memcpy (S->rgb_high, rgb, 4, double); S->z_high = z_cut; } /* Recompute differences used in gmt_get_rgb_from_z */ for (k = 0; k < 4; k++) S->rgb_diff[k] = S->rgb_high[k] - S->rgb_low[k]; for (k = 0; k < 4; k++) S->hsv_diff[k] = S->hsv_high[k] - S->hsv_low[k]; /* Recompute inverse stepsize */ S->i_dz = 1.0 / (S->z_high - S->z_low); } bool gmt_consider_current_cpt (struct GMTAPI_CTRL *API, bool *active, char **arg) { /* Modern mode only: Detect if no CPT is given but -C was set. * If -C[+u|U] is given (i.e., no cpt) then we take that to mean * we want to use the current CPT. If no current CPT exist then we do nothing * and the module with fail or create a CPT as needed. */ char *cpt = NULL; bool ret = true; if (API->GMT->current.setting.run_mode == GMT_CLASSIC) return false; /* No such thing in classic mode */ if (!(*active)) return false; /* Did not give that option so nothing to do */ if (arg == NULL) return false; /* No text pointer to work with */ if (gmt_M_cpt_mod (*arg)) { /* Gave modifiers for a unit change) */ char string[PATH_MAX] = {""}; if ((cpt = gmt_get_current_item (API->GMT, "cpt", false)) == NULL) return false; /* No current CPT */ snprintf (string, PATH_MAX, "%s%s", cpt, *arg); /* Append the modifiers to the current CPT name */ gmt_M_str_free (cpt); gmt_M_str_free (*arg); *arg = strdup (string); /* Pass back the name of the current CPT with modifiers */ } else if (*arg == NULL) { /* Noting given */ if ((cpt = gmt_get_current_item (API->GMT, "cpt", false)) == NULL) return false; /* No current CPT */ *arg = cpt; /* Pass back the name of the current CPT */ } else /* Got something already */ ret = false; return ret; } unsigned int gmt_set_interpolate_mode (struct GMT_CTRL *GMT, unsigned int mode, unsigned int type) { /* Convenience function that hides away the embedding of mode and type via the GMT_SPLINE_SLOPE factor */ gmt_M_unused (GMT); return (mode + GMT_SPLINE_SLOPE * type); } /*! . */ GMT_LOCAL double gmtsupport_csplint (struct GMT_CTRL *GMT, double *x, double *y, double *c, double xp, uint64_t klo) { uint64_t khi; double h, ih, b, a, yp; gmt_M_unused(GMT); khi = klo + 1; h = x[khi] - x[klo]; ih = 1.0 / h; a = (x[khi] - xp) * ih; b = (xp - x[klo]) * ih; yp = a * y[klo] + b * y[khi] + ((a*a*a - a) * c[klo] + (b*b*b - b) * c[khi]) * (h*h) / 6.0; return (yp); } /*! . */ GMT_LOCAL double gmtsupport_csplintslp (struct GMT_CTRL *GMT, double *x, double *y, double *c, double xp, uint64_t klo) { /* As gmtsupport_csplint but returns the first derivative instead */ uint64_t khi; double h, ih, b, a, dypdx; gmt_M_unused(GMT); khi = klo + 1; h = x[khi] - x[klo]; ih = 1.0 / h; a = (x[khi] - xp) * ih; b = (xp - x[klo]) * ih; dypdx = ih * (y[khi] - y[klo]) + ((3.0*b*b - 1.0) * c[khi] - (3.0*a*a - 1.0) * c[klo]) * h / 6.0; return (dypdx); } /*! . */ GMT_LOCAL double gmtsupport_csplintcurv (struct GMT_CTRL *GMT, double *x, double *c, double xp, uint64_t klo) { /* As gmtsupport_csplint but returns the 2nd derivative instead */ uint64_t khi; double h, ih, b, a, d2ypdx2; gmt_M_unused(GMT); khi = klo + 1; h = x[khi] - x[klo]; ih = 1.0 / h; a = (x[khi] - xp) * ih; b = (xp - x[klo]) * ih; d2ypdx2 = -(b * c[khi] + a * c[klo]); return (d2ypdx2); } /*! . */ int gmtlib_smooth_spline (struct GMT_CTRL *GMT, double *x, double *y, double *w, double p, uint64_t n, int deriv, double *a) { /* Implements a smoothing splines as per Lancaster & Salkauskas [1980]. * w is optional weights or NULL, p is fit parameter. The x values are monotonic */ int error; uint64_t n3 = n - 3, n2 = n - 2, n1 = n - 1, k, kk, kkt; double ip = 1.0 / p, sum_w = 0; double *D = NULL, *Dt = NULL, *B = NULL, *K = NULL, *Q = NULL; double *lambda = NULL, *dtau = NULL, *inv_dtau = NULL, *c = NULL, *s = NULL; /* Create the inverse lambda diagonal vector (assume w[k] = sigma_k) */ lambda = gmt_M_memory (GMT, NULL, n, double); for (k = 0; k < n; k++) { lambda[k] = (w == NULL) ? 1.0 : w[k] * w[k]; sum_w += lambda[k]; } for (k = 0; k < n; k++) { lambda[k] /= sum_w; /* Normalize weights to sum to 1 */ lambda[k] *= ip; } /* Create dtau and inv_dtau vectors */ dtau = gmt_M_memory (GMT, NULL, n1, double); inv_dtau = gmt_M_memory (GMT, NULL, n1, double); for (k = 0; k < n1; k++) { dtau[k] = x[k+1] - x[k]; inv_dtau[k] = 1.0 / dtau[k]; } /* Create B matrix */ B = gmt_M_memory (GMT, NULL, n2*n2, double); for (k = 0; k < n2; k++) { /* For rows 0 to n-3, a total of n-2 rows in B */ kk = k * n2 + k; /* Diagonal index */ B[kk] = (dtau[k] + dtau[k+1]) / 3.0; if (k > 0) B[kk-1] = dtau[k] / 6.0; /* All but first row has entry to left of diagonal */ if (k < n3) B[kk+1] = dtau[k+1] / 6.0; /* All but last row has entry to right of diagonal */ } /* Create D and Dt matrix and temp matrix K = inv_lambd * Dt */ D = gmt_M_memory (GMT, NULL, n2*n, double); Dt = gmt_M_memory (GMT, NULL, n*n2, double); K = gmt_M_memory (GMT, NULL, n2*n, double); for (k = 0; k < n2; k++) { /* For rows 0 to n-3, a total of n-2 rows in D */ kk = k * n + k; /* Diagonal index in D */ kkt = k * n2 + k; /* Diagonal index in transpose */ D[kk] = Dt[kkt] = inv_dtau[k]; D[kk+1] = Dt[kkt+n2] = -(inv_dtau[k] + inv_dtau[k+1]); D[kk+2] = Dt[kkt+2*n2] = inv_dtau[k+1]; K[kkt] = Dt[kkt] * lambda[k]; K[kkt+n2] = Dt[kkt+n2] * lambda[k+1]; K[kkt+2*n2] = Dt[kkt+2*n2] * lambda[k+2]; } /* Need to compute Q = D * K */ Q = gmt_M_memory (GMT, NULL, n2*n2, double); gmt_matrix_matrix_mult (GMT, D, K, n2, n, n2, Q); /* Now add B to Q */ gmt_matrix_matrix_add (GMT, Q, B, n2, n2, Q); /* Set c = D * y */ c = gmt_M_memory (GMT, NULL, n, double); gmt_matrix_matrix_mult (GMT, D, y, n2, n, 1U, c); /* Solve Q*x = c, with x returned in c; these are the curvatures (s" in the notes) */ error = gmt_gauss (GMT, Q, c, n2, n2, true); if (error) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Gauss returns error code %d\n", error); error = GMT_RUNTIME_ERROR; goto END_IT; } /* Set Q = K * c which is a vector of length n */ gmt_matrix_matrix_mult (GMT, K, c, n, n2, 1, Q); /* Solution for s = y - Q as in the notes */ s = gmt_M_memory (GMT, NULL, n, double); for (k = 0; k < n; k++) s[k] = y[k] - Q[k]; /* Adjust c (n-2 items) by adding the first and last c = conditions */ for (k = n2; k > 1; k--) c[k] = c[k-1]; c[0] = c[n-1] = 0.0; /* THe natural BCs */ /* Now build the coefficient matrix a (5.33) with y = s and c = s" */ for (k = kk = 0; k < n1; k++) { /* The internal segments */ a[kk++] = c[k] * inv_dtau[k] / 6.0; a[kk++] = c[k+1] * inv_dtau[k] / 6.0; a[kk++] = s[k] * inv_dtau[k] - c[k] * dtau[k] / 6.0; a[kk++] = s[k+1] * inv_dtau[k] - c[k+1] * dtau[k] / 6.0; } if (deriv == 1) { /* Take first derivative and update coefficients */ for (k = kk = 0; k < n1; k++) { /* The internal segments */ a[kk++] *= -3; a[kk++] *= +3; a[kk] = a[kk+1] - a[kk]; kk++; a[kk++] = 0.0; } } else if (deriv == 2) { /* Take second derivative and update coefficients */ for (k = kk = 0; k < n1; k++) { /* The internal segments */ a[kk++] *= 6; a[kk++] *= 6; a[kk++] = 0.0; a[kk++] = 0.0; } } END_IT: gmt_M_free (GMT, c); gmt_M_free (GMT, s); gmt_M_free (GMT, lambda); gmt_M_free (GMT, dtau); gmt_M_free (GMT, inv_dtau); gmt_M_free (GMT, Q); gmt_M_free (GMT, K); gmt_M_free (GMT, D); gmt_M_free (GMT, Dt); gmt_M_free (GMT, B); return (error); } /*! . */ GMT_LOCAL int gmtsupport_intpol_sub (struct GMT_CTRL *GMT, double *x, double *y, double *w, uint64_t n, uint64_t m, double *u, double *v, double p, int mode) { /* Does the main work of interpolating a section that has no NaNs */ uint64_t i, j; int err_flag = 0, spline_mode = mode % GMT_SPLINE_SLOPE, deriv = mode / GMT_SPLINE_SLOPE; double dx, dx1, dx2, x_min, x_max, *c = NULL; /* Set minimum and maximum */ if (spline_mode == GMT_SPLINE_NN) { x_min = (3*x[0] - x[1]) / 2; x_max = (3*x[n-1] - x[n-2]) / 2; } else { x_min = x[0]; x_max = x[n-1]; } /* Allocate memory for spline factors */ if (spline_mode == GMT_SPLINE_AKIMA) { /* Akima's spline */ c = gmt_M_memory (GMT, NULL, 3*n, double); err_flag = gmtlib_akima (GMT, x, y, n, c); } else if (spline_mode == GMT_SPLINE_CUBIC) { /* Natural cubic spline */ c = gmt_M_memory (GMT, NULL, 3*n, double); err_flag = gmtlib_cspline (GMT, x, y, n, c); } else if (spline_mode == GMT_SPLINE_SMOOTH) { /* Smoothing spline */ c = gmt_M_memory (GMT, NULL, 4*n, double); err_flag = gmtlib_smooth_spline (GMT, x, y, w, p, n, deriv, c); } if (err_flag != 0) { gmt_M_free (GMT, c); return (err_flag); } /* Compute the interpolated values from the coefficients */ j = 0; for (i = 0; i < m; i++) { if (u[i] < x_min || u[i] > x_max) { /* Desired point outside data range */ if (GMT->current.setting.extrapolate_val[0] == GMT_EXTRAPOLATE_NONE) { v[i] = GMT->session.d_NaN; continue; } else if (GMT->current.setting.extrapolate_val[0] == GMT_EXTRAPOLATE_CONSTANT) { v[i] = GMT->current.setting.extrapolate_val[1]; continue; } } while (j > 0 && x[j] > u[i]) j--; /* In case u is not sorted */ while (j < n && x[j] <= u[i]) j++; if (j == n) j--; if (j > 0) j--; switch (mode) { case GMT_SPLINE_LINEAR: /* Linear spline v(u) */ dx = u[i] - x[j]; v[i] = (y[j+1]-y[j])*dx/(x[j+1]-x[j]) + y[j]; break; case GMT_SPLINE_AKIMA: /* Akime's spline u(v) */ dx = u[i] - x[j]; v[i] = ((c[3*j+2]*dx + c[3*j+1])*dx + c[3*j])*dx + y[j]; break; case GMT_SPLINE_CUBIC: /* Natural cubic spline u(v) */ v[i] = gmtsupport_csplint (GMT, x, y, c, u[i], j); break; case GMT_SPLINE_NN: /* Nearest neighbor value */ v[i] = ((u[i] - x[j]) < (x[j+1] - u[i])) ? y[j] : y[j+1]; break; case GMT_SPLINE_LINEAR+GMT_SPLINE_SLOPE: /* Linear spline v'(u) */ v[i] = (y[j+1]-y[j])/(x[j+1]-x[j]); break; case GMT_SPLINE_AKIMA+GMT_SPLINE_SLOPE: /* Akime's spline u'(v) */ dx = u[i] - x[j]; v[i] = (3.0*c[3*j+2]*dx + 2.0*c[3*j+1])*dx + c[3*j]; break; case GMT_SPLINE_CUBIC+GMT_SPLINE_SLOPE: /* Natural cubic spline u'(v) */ v[i] = gmtsupport_csplintslp (GMT, x, y, c, u[i], j); break; case GMT_SPLINE_NN+GMT_SPLINE_SLOPE: /* Nearest neighbor slopes are zero */ v[i] = 0.0; break; case GMT_SPLINE_LINEAR+GMT_SPLINE_CURVATURE: /* Linear spline v"(u) */ v[i] = 0.0; break; case GMT_SPLINE_AKIMA+GMT_SPLINE_CURVATURE: /* Akime's spline u"(v) */ dx = u[i] - x[j]; v[i] = 6.0*c[3*j+2]*dx + 2.0*c[3*j+1]; break; case GMT_SPLINE_CUBIC+GMT_SPLINE_CURVATURE: /* Natural cubic spline u"(v) */ v[i] = gmtsupport_csplintcurv (GMT, x, c, u[i], j); break; case GMT_SPLINE_NN+GMT_SPLINE_CURVATURE: /* Nearest neighbor curvatures are zero */ v[i] = 0.0; break; case GMT_SPLINE_SMOOTH: /* Smoothing spline */ dx1 = x[j+1] - u[i]; dx2 = u[i] - x[j]; v[i] = c[4*j] * pow (dx1, 3.0) + c[4*j+1] * pow (dx2, 3.0) + c[4*j+2] * dx1 + c[4*j+3] * dx2; break; case GMT_SPLINE_SMOOTH+GMT_SPLINE_SLOPE: /* Derivative of smoothing spline */ dx1 = x[j+1] - u[i]; dx2 = u[i] - x[j]; v[i] = c[4*j] * pow (dx1, 2.0) + c[4*j+1] * pow (dx2, 2.0) + c[4*j+2]; break; case GMT_SPLINE_SMOOTH+GMT_SPLINE_CURVATURE: /* Curvature of smoothing spline */ dx1 = x[j+1] - u[i]; dx2 = u[i] - x[j]; v[i] = c[4*j] * dx1 + c[4*j+1] * dx2; break; } } gmt_M_free (GMT, c); return (GMT_NOERROR); } /*! . */ GMT_LOCAL void gmtsupport_intpol_reverse (double *x, double *u, uint64_t n, uint64_t m) { /* Changes sign on x and u */ uint64_t i; for (i = 0; i < n; i++) x[i] = -x[i]; for (i = 0; i < m; i++) u[i] = -u[i]; } /*! Non-square matrix transpose of matrix of size r x c and base address A */ GMT_LOCAL unsigned int gmtsupport_is_set (unsigned int *mark, uint64_t k, unsigned int *bits) { /* Return nonzero if this bit is set */ uint64_t w = k / 32ULL; unsigned int b = (unsigned int)(k % 32ULL); return (mark[w] & bits[b]); } /*! . */ GMT_LOCAL void gmtsupport_set_bit (unsigned int *mark, uint64_t k, unsigned int *bits) { /* Set this bit */ uint64_t w = k / 32ULL; unsigned int b = (unsigned int)(k % 32ULL); mark[w] |= bits[b]; } /*! . */ GMT_LOCAL void gmtsupport_decorate_free (struct GMT_CTRL *GMT, struct GMT_DECORATE *G) { /* Free memory used by decorate */ GMT_Destroy_Data (GMT->parent, &(G->X)); if (G->f_n) { /* Array for fixed points */ gmt_M_free (GMT, G->f_xy[GMT_X]); gmt_M_free (GMT, G->f_xy[GMT_Y]); } } /*! . */ GMT_LOCAL int gmtsupport_code_to_lonlat (struct GMT_CTRL *GMT, char *code, double *lon, double *lat) { int i, n, error = 0; bool z_OK = false; n = (int)strlen (code); if (n != 2) return (1); for (i = 0; i < 2; i++) { switch (code[i]) { case 'l': case 'L': /* Left */ *lon = GMT->common.R.wesn[XLO]; break; case 'c': case 'C': /* center */ *lon = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]); break; case 'r': case 'R': /* right */ *lon = GMT->common.R.wesn[XHI]; break; case 'b': case 'B': /* bottom */ *lat = GMT->common.R.wesn[YLO]; break; case 'm': case 'M': /* center */ *lat = 0.5 * (GMT->common.R.wesn[YLO] + GMT->common.R.wesn[YHI]); break; case 't': case 'T': /* top */ *lat = GMT->common.R.wesn[YHI]; break; case 'z': case 'Z': /* z-value */ z_OK = true; break; case '+': /* zmax-location */ if (z_OK) *lon = *lat = DBL_MAX; else error++; break; case '-': /* zmin-location */ if (z_OK) *lon = *lat = -DBL_MAX; else error++; break; default: error++; break; } } if (error) GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized location code %s\n", code); return (error); } /*! . */ GMT_LOCAL double gmtsupport_determine_endpoint (struct GMT_CTRL *GMT, double x0, double y0, double length, double az, double *x1, double *y1) { /* compute point a distance length from origin along azimuth, return point separation */ double s_az, c_az; sincosd (az, &s_az, &c_az); if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Spherical solution */ double s0, c0, s_d, c_d; double d = (length / GMT->current.map.dist[GMT_MAP_DIST].scale); /* Convert from chosen unit to meter or degree */ if (!GMT->current.map.dist[GMT_MAP_DIST].arc) d /= GMT->current.proj.DIST_M_PR_DEG; /* Convert meter to spherical degrees */ sincosd (y0, &s0, &c0); sincosd (d, &s_d, &c_d); *x1 = x0 + atand (s_d * s_az / (c0 * c_d - s0 * s_d * c_az)); *y1 = d_asind (s0 * c_d + c0 * s_d * c_az); } else { /* Cartesian solution */ *x1 = x0 + length * s_az; *y1 = y0 + length * c_az; } return (gmt_distance (GMT, x0, y0, *x1, *y1)); } /*! . */ GMT_LOCAL double gmtsupport_determine_endpoints (struct GMT_CTRL *GMT, double x[], double y[], double length, double az) { double s_az, c_az; sincosd (az, &s_az, &c_az); length /= 2.0; /* Going half-way in each direction */ if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Spherical solution */ double s0, c0, s_d, c_d; double d = (length / GMT->current.map.dist[GMT_MAP_DIST].scale); /* Convert from chosen unit to meter or degree */ if (!GMT->current.map.dist[GMT_MAP_DIST].arc) d /= GMT->current.proj.DIST_M_PR_DEG; /* Convert meter to spherical degrees */ sincosd (y[0], &s0, &c0); sincosd (d, &s_d, &c_d); x[1] = x[0] + atand (s_d * s_az / (c0 * c_d - s0 * s_d * c_az)); y[1] = d_asind (s0 * c_d + c0 * s_d * c_az); /* Then redo start point by going in the opposite direction */ sincosd (az+180.0, &s_az, &c_az); x[0] = x[0] + atand (s_d * s_az / (c0 * c_d - s0 * s_d * c_az)); y[0] = d_asind (s0 * c_d + c0 * s_d * c_az); } else { /* Cartesian solution */ x[1] = x[0] + length * s_az; y[1] = y[0] + length * c_az; x[0] = x[0] - length * s_az; y[0] = y[0] - length * c_az; } return (gmt_distance (GMT, x[0], y[0], x[1], y[1])); } /*! . */ GMT_LOCAL uint64_t gmtsupport_determine_circle (struct GMT_CTRL *GMT, double x0, double y0, double r, double x[], double y[], unsigned int n) { /* Given an origin, radius, and n points, compute a circular path and return it */ unsigned int k; uint64_t np = n; double d_angle = 360.0 / n; if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Spherical solution */ double R[3][3], v[3], v_prime[3], lat; double colat = (r / GMT->current.map.dist[GMT_MAP_DIST].scale); /* Convert from chosen unit to meter or degree */ if (!GMT->current.map.dist[GMT_MAP_DIST].arc) colat /= GMT->current.proj.DIST_M_PR_DEG; /* Convert meter to spherical degrees */ lat = 90.0 - colat; /* Get small circle latitude about NP */ /* Set up small circle in local coordinates, convert to Cartesian, rotate, and covert to geo */ /* Rotation pole is 90 away from x0, at equator, and rotates by 90-y0 degrees */ gmt_make_rot_matrix (GMT, x0 + 90.0, 0.0, 90.0 - y0, R); for (k = 0; k < n; k++) { gmt_geo_to_cart (GMT, lat, d_angle * k, v, true); gmt_matrix_vect_mult (GMT, 3U, R, v, v_prime); gmt_cart_to_geo (GMT, &y[k], &x[k], v_prime, true); } } else { /* Cartesian solution */ double s_az, c_az; for (k = 0; k < n; k++) { sincosd (d_angle * k, &s_az, &c_az); x[k] = x0 + r * s_az; y[k] = y0 + r * c_az; } } return (np); } /*! . */ GMT_LOCAL void gmtsupport_line_angle_ave (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, int half, int angle_type, struct GMT_LABEL *L) { int64_t j, sstart, sstop, nn; double sum_x2 = 0.0, sum_xy = 0.0, sum_y2 = 0.0, dx, dy; if (start == stop) { /* Can happen if we want no smoothing but landed exactly on a knot point */ if (start > 0) start--; else if (stop < (n-1)) stop++; } sstart = start - half; sstop = stop + half; nn = n; for (j = sstart; j <= sstop; j++) { /* L2 fit for slope over this range of points */ if (j < 0 || j >= nn) continue; dx = x[j] - L->x; dy = y[j] - L->y; sum_x2 += dx * dx; sum_y2 += dy * dy; sum_xy += dx * dy; } if (sum_y2 < GMT_CONV8_LIMIT) /* Line is horizontal */ L->line_angle = 0.0; else if (sum_x2 < GMT_CONV8_LIMIT) /* Line is vertical */ L->line_angle = 90.0; else L->line_angle = (gmt_M_is_zero (sum_xy)) ? 90.0 : d_atan2d (sum_xy, sum_x2); if (angle_type == 2) { /* Just return the fixed angle given (unless NaN) */ if (gmt_M_is_dnan (cangle)) /* Cannot use this angle - default to along-line angle */ angle_type = 0; else L->angle = cangle; } if (angle_type != 2) { /* Must base label angle on the contour angle */ L->angle = L->line_angle + angle_type * 90.0; /* May add 90 to get normal */ if (L->angle < 0.0) L->angle += 360.0; if (L->angle > 90.0 && L->angle < 270) L->angle -= 180.0; } GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Ave: Label Line angle = %g start/stop = %d/%d atan2d (%g, %g) Label angle = %g\n", L->line_angle, (int)start, (int)stop, sum_xy, sum_x2, L->angle); } /*! . */ GMT_LOCAL void gmtsupport_line_angle_line (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, int angle_type, struct GMT_LABEL *L) { double dx, dy; if (start == stop) { /* Can happen if we want no smoothing but landed exactly on a knot point */ if (start > 0) start--; else if (stop < (n-1)) stop++; } if (stop >= n) stop = n-1; dx = x[stop] - x[start]; dy = y[stop] - y[start]; L->line_angle = d_atan2d (dy, dx); if (angle_type == 2) { /* Just return the fixed angle given (unless NaN) */ if (gmt_M_is_dnan (cangle)) /* Cannot use this angle - default to along-line angle */ angle_type = 0; else L->angle = cangle; } if (angle_type != 2) { /* Must base label angle on the contour angle */ L->angle = L->line_angle + angle_type * 90.0; /* May add 90 to get normal */ if (L->angle < 0.0) L->angle += 360.0; if (L->angle > 90.0 && L->angle < 270) L->angle -= 180.0; } GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Vec: Label Line angle = %g start/stop = %d/%d atan2d (%g, %g) Label angle = %g\n", L->line_angle, (int)start, (int)stop, dy, dx, L->angle); } /*! . */ GMT_LOCAL void gmtsupport_contlabel_angle_ave (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, struct GMT_LABEL *L, struct GMT_CONTOUR *G) { gmtsupport_line_angle_ave (GMT, x, y, start, stop, cangle, n, G->half_width, G->angle_type, L); } /*! . */ GMT_LOCAL void gmtsupport_contlabel_angle_line (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, struct GMT_LABEL *L, struct GMT_CONTOUR *G) { gmtsupport_line_angle_line (GMT, x, y, start, stop, cangle, n, G->angle_type, L); } /*! . */ GMT_LOCAL void gmtsupport_contlabel_angle (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, bool contour, struct GMT_LABEL *L, struct GMT_CONTOUR *G) { /* Sets L->line_angle and L->angle */ if ((G->nudge_flag == 2 && G->half_width == UINT_MAX ) || G->half_width == 0) { /* Want line-angle to follow line */ gmtsupport_contlabel_angle_line (GMT, x, y, start, stop, cangle, n, L, G); } else if (G->half_width == UINT_MAX) { /* Automatic width specification */ /* Try to come up with a number that is small for short lines and grow slowly for larger lines */ G->half_width = MAX (1, lrint (ceil (log10 (0.3333333333*n)))); G->half_width *= G->half_width; /* New guess at half-width */ GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Automatic label-averaging half_width = %d [n = %d]\n", G->half_width, (int)n); if (G->half_width == 1) gmtsupport_contlabel_angle_line (GMT, x, y, start, stop, cangle, n, L, G); else gmtsupport_contlabel_angle_ave (GMT, x, y, start, stop, cangle, n, L, G); G->half_width = UINT_MAX; /* Reset back to auto */ } else { /* Go width the selected half-width */ gmtsupport_contlabel_angle_ave (GMT, x, y, start, stop, cangle, n, L, G); } if (contour) { /* Limit line_angle to -90/+90 */ if (L->line_angle > +90.0) L->line_angle -= 180.0; if (L->line_angle < -90.0) L->line_angle += 180.0; } } /*! . */ GMT_LOCAL void gmtsupport_decorated_angle_ave (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, struct GMT_LABEL *L, struct GMT_DECORATE *G) { gmtsupport_line_angle_ave (GMT, x, y, start, stop, cangle, n, G->half_width, G->angle_type, L); } /*! . */ GMT_LOCAL void gmtsupport_decorated_angle_line (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, struct GMT_LABEL *L, struct GMT_DECORATE *G) { gmtsupport_line_angle_line (GMT, x, y, start, stop, cangle, n, G->angle_type, L); } /*! . */ GMT_LOCAL void gmtsupport_decorated_angle (struct GMT_CTRL *GMT, double x[], double y[], uint64_t start, uint64_t stop, double cangle, uint64_t n, bool contour, struct GMT_LABEL *L, struct GMT_DECORATE *G) { /* Sets L->line_angle and L->angle */ if ((G->nudge_flag == 2 && G->half_width == UINT_MAX ) || G->half_width == 0) { /* Want line-angle to follow line */ gmtsupport_decorated_angle_line (GMT, x, y, start, stop, cangle, n, L, G); } else if (G->half_width == UINT_MAX) { /* Automatic width specification */ /* Try to come up with a number that is small for short lines and grow slowly for larger lines */ G->half_width = MAX (1, lrint (ceil (log10 (0.3333333333*n)))); G->half_width *= G->half_width; /* New guess at half-width */ GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Automatic label-averaging half_width = %d [n = %d]\n", G->half_width, (int)n); if (G->half_width == 1) gmtsupport_decorated_angle_line (GMT, x, y, start, stop, cangle, n, L, G); else gmtsupport_decorated_angle_ave (GMT, x, y, start, stop, cangle, n, L, G); G->half_width = UINT_MAX; /* Reset back to auto */ } else { /* Go width the selected half-width */ gmtsupport_decorated_angle_ave (GMT, x, y, start, stop, cangle, n, L, G); } if (contour) { /* Limit line_angle to -90/+90 */ if (L->line_angle > +90.0) L->line_angle -= 180.0; if (L->line_angle < -90.0) L->line_angle += 180.0; } } /*! . */ GMT_LOCAL int gmtsupport_sort_label_struct (const void *p_1, const void *p_2) { const struct GMT_LABEL **point_1 = (const struct GMT_LABEL **)p_1, **point_2 = (const struct GMT_LABEL **)p_2; if ((*point_1)->dist < (*point_2)->dist) return -1; if ((*point_1)->dist > (*point_2)->dist) return +1; return 0; } /*! . */ GMT_LOCAL void gmtsupport_contlabel_fixpath (struct GMT_CTRL *GMT, double **xin, double **yin, double d[], uint64_t *n, struct GMT_CONTOUR *G) { /* Sorts labels based on distance and inserts the label (x,y) point into the x,y path */ uint64_t i, j, k, np; double *xp = NULL, *yp = NULL, *x = NULL, *y = NULL; if (G->n_label == 0) return; /* No labels, no need to insert points */ /* Sort labels based on distance along contour if more than 1 */ if (G->n_label > 1) qsort (G->L, G->n_label, sizeof (struct GMT_LABEL *), gmtsupport_sort_label_struct); np = *n + G->n_label; /* Length of extended path that includes inserted label coordinates */ xp = gmt_M_memory (GMT, NULL, np, double); yp = gmt_M_memory (GMT, NULL, np, double); x = *xin; y = *yin; /* Input coordinate arrays */ /* Make sure the xp, yp array contains the exact points at which labels should be placed */ for (k = 0, i = j = 0; i < *n && j < np && k < G->n_label; k++) { while (i < *n && d[i] < G->L[k]->dist) { /* Not at the label point yet - just copy points */ xp[j] = x[i]; yp[j] = y[i]; j++; i++; } /* Add the label point */ G->L[k]->node = j; /* Update node since we have been adding new points */ xp[j] = G->L[k]->x; yp[j] = G->L[k]->y; j++; } while (i < *n) { /* Append the rest of the path */ xp[j] = x[i]; yp[j] = y[i]; j++; i++; } gmt_M_free (GMT, x); /* Get rid of old path... */ gmt_M_free (GMT, y); *xin = xp; /* .. and pass out new path */ *yin = yp; *n = np; /* and the new length */ } /*! . */ GMT_LOCAL void gmtsupport_contlabel_addpath (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double zval, char *label, bool annot, struct GMT_CONTOUR *G) { uint64_t i; double s = 0.0, c = 1.0, sign = 1.0; struct GMT_CONTOUR_LINE *L = NULL; /* Adds this segment to the list of contour lines */ if (G->n_alloc == 0 || G->n_segments == G->n_alloc) { G->n_alloc = (G->n_alloc == 0) ? GMT_CHUNK : (G->n_alloc << 1); G->segment = gmt_M_memory (GMT, G->segment, G->n_alloc, struct GMT_CONTOUR_LINE *); } G->segment[G->n_segments] = gmt_M_memory (GMT, NULL, 1, struct GMT_CONTOUR_LINE); L = G->segment[G->n_segments]; /* Pointer to current segment */ L->n = n; L->x = gmt_M_memory (GMT, NULL, L->n, double); L->y = gmt_M_memory (GMT, NULL, L->n, double); gmt_M_memcpy (L->x, x, L->n, double); gmt_M_memcpy (L->y, y, L->n, double); gmt_M_memcpy (&L->pen, &G->line_pen, 1, struct GMT_PEN); /* gmt_M_memcpy (L->rgb,&G->rgb, 4, double); */ gmt_M_memcpy (L->rgb, &G->font_label.fill.rgb, 4, double); L->name = gmt_M_memory (GMT, NULL, strlen (label)+1, char); strcpy (L->name, label); L->annot = annot; L->z = zval; if (G->n_label) { /* There are labels */ L->n_labels = G->n_label; L->L = gmt_M_memory (GMT, NULL, L->n_labels, struct GMT_LABEL); for (i = 0; i < L->n_labels; i++) { L->L[i].x = G->L[i]->x; L->L[i].y = G->L[i]->y; L->L[i].line_angle = G->L[i]->line_angle; if (G->nudge_flag) { /* Must adjust point a bit */ if (G->nudge_flag == 2) sincosd (L->L[i].line_angle, &s, &c); /* If N+1 or N-1 is used we want positive x nudge to extend away from end point */ sign = (G->number_placement) ? (double)L->L[i].end : 1.0; L->L[i].x += sign * (G->nudge[GMT_X] * c - G->nudge[GMT_Y] * s); L->L[i].y += sign * (G->nudge[GMT_X] * s + G->nudge[GMT_Y] * c); } L->L[i].angle = G->L[i]->angle; L->L[i].dist = G->L[i]->dist; L->L[i].node = G->L[i]->node; L->L[i].label = gmt_M_memory (GMT, NULL, strlen (G->L[i]->label)+1, char); gmt_M_memcpy (L->L[i].rgb, G->L[i]->rgb, 4, double); strcpy (L->L[i].label, G->L[i]->label); } } G->n_segments++; } /*! . */ GMT_LOCAL void gmtsupport_get_radii_of_curvature (double x[], double y[], uint64_t n, double r[]) { /* Calculates radius of curvature along the spatial curve x(t), y(t) */ uint64_t i, im, ip; double a, b, c, d, e, f, x0, y0, denom; for (im = 0, i = 1, ip = 2; ip < n; i++, im++, ip++) { a = (x[im] - x[i]); b = (y[im] - y[i]); e = 0.5 * (x[im] * x[im] + y[im] * y[im] - x[i] * x[i] - y[i] * y[i]); c = (x[i] - x[ip]); d = (y[i] - y[ip]); f = 0.5 * (x[i] * x[i] + y[i] * y[i] - x[ip] * x[ip] - y[ip] * y[ip]); denom = b * c - a * d; if (denom == 0.0) r[i] = DBL_MAX; else { x0 = (-d * e + b * f) / denom; y0 = (c * e - a * f) / denom; r[i] = hypot (x[i] - x0, y[i] - y0); } } r[0] = r[n-1] = DBL_MAX; /* Boundary conditions has zero curvature at end points so r = inf */ } /*! . */ GMT_LOCAL void gmtsupport_edge_contour (struct GMT_CTRL *GMT, struct GMT_GRID *G, unsigned int col, unsigned int row, unsigned int side, double d, double *x, double *y) { gmt_M_unused(GMT); if (side == 0) { *x = gmt_M_grd_col_to_x (GMT, col+d, G->header); *y = gmt_M_grd_row_to_y (GMT, row, G->header); } else if (side == 1) { *x = gmt_M_grd_col_to_x (GMT, col+1, G->header); *y = gmt_M_grd_row_to_y (GMT, row-d, G->header); } else if (side == 2) { *x = gmt_M_grd_col_to_x (GMT, col+1-d, G->header); *y = gmt_M_grd_row_to_y (GMT, row-1, G->header); } else { *x = gmt_M_grd_col_to_x (GMT, col, G->header); *y = gmt_M_grd_row_to_y (GMT, row-1+d, G->header); } } /*! . */ GMT_LOCAL void gmtsupport_setcontjump (gmt_grdfloat *z, uint64_t nz) { /* This routine will check if there is a 360 jump problem * among these coordinates and adjust them accordingly so * that subsequent testing can determine if a zero contour * goes through these edges. E.g., values like 359, 1 * should become -1, 1 after this function */ uint64_t i; bool jump = false; gmt_grdfloat dz; #ifdef DOUBLE_PRECISION_GRID for (i = 1; !jump && i < nz; i++) { dz = z[i] - z[0]; if (fabs (dz) > 180.0) jump = true; } if (!jump) return; z[0] = fmod (z[0], 360.0); if (z[0] > 180.0) z[0] -= 360.0; else if (z[0] < -180.0) z[0] += 360.0; for (i = 1; i < nz; i++) { if (z[i] > 180.0) z[i] -= 360.0; else if (z[i] < -180.0) z[i] += 360.0; dz = z[i] - z[0]; if (fabs (dz) > 180.0) z[i] -= copysign (360.0, dz); z[i] = fmod (z[i], 360.0); } #else for (i = 1; !jump && i < nz; i++) { dz = z[i] - z[0]; if (fabsf (dz) > 180.0f) jump = true; } if (!jump) return; z[0] = fmodf (z[0], 360.0f); if (z[0] > 180.0f) z[0] -= 360.0f; else if (z[0] < -180.0f) z[0] += 360.0f; for (i = 1; i < nz; i++) { if (z[i] > 180.0f) z[i] -= 360.0f; else if (z[i] < -180.0f) z[i] += 360.0f; dz = z[i] - z[0]; if (fabsf (dz) > 180.0f) z[i] -= copysignf (360.0f, dz); z[i] = fmodf (z[i], 360.0f); } #endif } /*! . */ GMT_LOCAL uint64_t gmtsupport_trace_contour (struct GMT_CTRL *GMT, struct GMT_GRID *G, bool test, unsigned int *edge, double **x, double **y, unsigned int col, unsigned int row, unsigned int side, uint64_t offset, unsigned int *bit, unsigned int *nan_flag) { /* Note: side must be signed due to calculations like (side-2)%2 which will not work with unsigned */ unsigned int side_in, this_side, old_side, n_exits, opposite_side, n_nan, edge_word, edge_bit, periodic = 0; int p[5] = {0, 0, 0, 0, 0}, mx; bool more, crossed = false; size_t n_alloc; uint64_t n = 1, m, ij0, ij_in, ij; gmt_grdfloat z[5] = {0.0, 0.0, 0.0, 0.0, 0.0}, dz; double xk[5] = {0.0, 0.0, 0.0, 0.0, 0.0}, dist1, dist2, *xx = NULL, *yy = NULL; static int d_col[5] = {0, 1, 0, 0, 0}, d_row[5] = {0, 0, -1, 0, 0}, d_side[5] = {0, 1, 0, 1, 0}; #if 0 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (G->header); /* WORK IN PROGRESS. Currently has trouble with some contours streaking across a W/E map */ /* For periodic grids we must also be sensitive to exiting in west and reenter in east, and vice versa */ if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT) periodic = 1; else if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_NOREPEAT) periodic = 2; #endif /* Note: We need gmt_M_ijp to get us the index into the padded G->data whereas we use gmt_M_ij0 to get the corresponding index for non-padded edge array */ mx = G->header->mx; /* Need signed mx below */ p[0] = p[4] = 0; p[1] = 1; p[2] = 1 - mx; p[3] = -mx; /* Padded offsets for G->data array */ *nan_flag = 0; /* Check if this edge was already used */ ij0 = gmt_M_ij0 (G->header, row + d_row[side], col + d_col[side]); /* Index to use with the edge array */ edge_word = (unsigned int)(ij0 / 32 + d_side[side] * offset); edge_bit = (unsigned int)(ij0 % 32); if (test && (edge[edge_word] & bit[edge_bit])) return (0); /* Been here */ ij_in = gmt_M_ijp (G->header, row, col); /* Index to use with the padded G->data array */ side_in = side; /* First check if contour cuts the starting edge */ z[0] = G->data[ij_in+p[side]]; z[1] = G->data[ij_in+p[side+1]]; if (GMT->current.map.z_periodic) gmtsupport_setcontjump (z, 2); if (!(z[0] * z[1] < 0.0f)) return (0); /* This formulation will also return if one of the z's is NaN */ n_alloc = GMT_CHUNK; m = n_alloc - 2; xx = gmt_M_memory (GMT, NULL, n_alloc, double); yy = gmt_M_memory (GMT, NULL, n_alloc, double); gmtsupport_edge_contour (GMT, G, col, row, side, z[0] / (z[0] - z[1]), &(xx[0]), &(yy[0])); edge[edge_word] |= bit[edge_bit]; more = true; do { ij = gmt_M_ijp (G->header, row, col); /* Index to use with the padded G->data array */ /* If this is the box and edge we started from, explicitly close the polygon and exit */ if (n > 1 && ij == ij_in && side == side_in) { /* Yes, we close and exit */ xx[n-1] = xx[0]; yy[n-1] = yy[0]; /* Ensure no roundoff */ more = false; continue; } n_exits = 0; old_side = side; for (this_side = 0; this_side < 5; this_side++) z[this_side] = G->data[ij+p[this_side]]; /* Copy the 4 corners to the z array and duplicate the 1st as a '5th' corner */ if (GMT->current.map.z_periodic) gmtsupport_setcontjump (z, 5); for (this_side = n_nan = 0; this_side < 4; this_side++) { /* Loop over the 4 box sides and count possible exits */ /* Skip if NaN encountered */ if (gmt_M_is_fnan (z[this_side+1]) || gmt_M_is_fnan (z[this_side])) { /* Check each side, defined by two nodes, for NaNs */ n_nan++; continue; } /* Skip if no zero-crossing on this edge */ if (z[this_side+1] * z[this_side] > 0.0f) continue; if ((dz = z[this_side] - z[this_side+1]) == 0.0f) continue; /* Save normalized distance along edge from corner this_side to crossing of edge this_side */ xk[this_side] = z[this_side] / dz; /* Skip if this is the entry edge of the current box (this_side == old_side) */ if (this_side == old_side) continue; /* Here we have a new crossing */ side = this_side; n_exits++; } xk[4] = xk[0]; /* Repeat the first value */ if (n > m) { /* Must try to allocate more memory */ n_alloc <<= 1; m = n_alloc - 2; xx = gmt_M_memory (GMT, xx, n_alloc, double); yy = gmt_M_memory (GMT, yy, n_alloc, double); } /* These are the possible outcomes of n_exits: 0: No exits. We must have struck a wall of NaNs 1: One exit. Take it! 2: Two exits is not possible! 3: Three exits means we have entered on a saddle point */ if (n_exits == 0) { /* We have hit a field of NaNs, finish contour */ more = false; *nan_flag = n_nan; continue; } else if (n_exits == 3) { /* Saddle point: Turn to the correct side */ opposite_side = (old_side + 2) % 4; /* Opposite side */ dist1 = xk[old_side] + xk[opposite_side]; dist2 = xk[old_side+1] + xk[opposite_side+1]; if (dist1 == dist2) { /* Perfect saddle, must use fixed saddle decision to avoid a later crossing */ /* Connect S with W and E with N in this case */ if (old_side == 0) side = 3; else if (old_side == 1) side = 2; else if (old_side == 2) side = 1; else side = 0; } else if (dist1 > dist2) side = (old_side + 1) % 4; else side = (old_side + 3) % 4; } gmtsupport_edge_contour (GMT, G, col, row, side, xk[side], &(xx[n]), &(yy[n])); n++; /* Mark the new edge as used */ ij0 = gmt_M_ij0 (G->header, row + d_row[side], col + d_col[side]); /* Non-padded index used for edge array */ edge_word = (unsigned int)(ij0 / 32 + d_side[side] * offset); edge_bit = (unsigned int)(ij0 % 32); edge[edge_word] |= bit[edge_bit]; if ((side == 0 && row == G->header->n_rows - 1) || (side == 2 && row == 1)) { /* Going out of grid at N or S */ more = false; continue; } else if ((side == 1 && col == G->header->n_columns - 2) || (side == 3 && col == 0)) { /* Going out of grid E-W */ if (periodic == 0) { /* And that is it. */ more = false; continue; } if (side == 1) { /* Exiting on east to reappear on west */ side = 3; /* So entering first col bin from left (3) */ col = 0; /* First column bin */ xx[n] = xx[n-1] - 360.0; } else { /* Exiting on west to reappear on east */ side = 1; /* So entering last col bin from right (1) */ col = G->header->n_columns - 2; /* Last column bin since next col is the boundary of the grid */ xx[n] = xx[n-1] + 360.0; } yy[n] = yy[n-1]; crossed = true; /* Flag that we crossed periodic boundary and add the duplicate coordinate point to array */ n++; ij0 = gmt_M_ij0 (G->header, row + d_row[side], col + d_col[side]); /* Flag this edge as used as well */ edge_word = (unsigned int)(ij0 / 32 + d_side[side] * offset); edge_bit = (unsigned int)(ij0 % 32); if (test && (edge[edge_word] & bit[edge_bit])) more = false; /* Been here before */ edge[edge_word] |= bit[edge_bit]; continue; } switch (side) { /* Move on to next box (col,row,side) */ case 0: row++; side = 2; break; /* Go to row below */ case 1: col++; side = 3; break; /* Go to col on right */ case 2: row--; side = 0; break; /* Go to row above */ case 3: col--; side = 1; break; /* Go to col on left */ } } while (more); xx = gmt_M_memory (GMT, xx, n, double); yy = gmt_M_memory (GMT, yy, n, double); if (crossed) { if ((gmtlib_determine_pole (GMT, xx, yy, n) == 0)) gmt_eliminate_lon_jumps (GMT, xx, n); /* Ensure longitudes are in the same quadrants */ } *x = xx; *y = yy; return (n); } /*! . */ GMT_LOCAL int64_t gmtsupport_smooth_contour (struct GMT_CTRL *GMT, double **x_in, double **y_in, uint64_t n, int sfactor, int stype) { /* Input n (x_in, y_in) points */ /* n_out = sfactor * n -1 */ /* Interpolation scheme used (0 = linear, 1 = Akima, 2 = Cubic spline, 3 = None */ uint64_t i, j, k, n_out; double ds, t_next, *x = NULL, *y = NULL; double *t_in = NULL, *t_out = NULL, *x_tmp = NULL, *y_tmp = NULL, x0, x1, y0, y1; bool *flag = NULL; if (sfactor == 0 || n < 4) return (n); /* Need at least 4 points to call Akima */ x = *x_in; y = *y_in; n_out = sfactor * n - 1; /* Number of new points */ t_in = gmt_M_memory (GMT, NULL, n, double); /* Create dummy distance values for input points, and throw out duplicate points while at it */ t_in[0] = 0.0; for (i = j = 1; i < n; i++) { if (gmt_M_is_geographic (GMT, GMT_IN) && gmt_M_360_range (x[i-1], x[i])) { ds = 0.0; /* 360 degree jumps are excluded */ } else ds = hypot ((x[i]-x[i-1]), (y[i]-y[i-1])); if (ds > 0.0) { t_in[j] = t_in[j-1] + ds; x[j] = x[i]; y[j] = y[i]; j++; } } n = j; /* May have lost some duplicates */ if (sfactor == 0 || n < 4) { /* Need at least 4 points to call Akima */ gmt_M_free (GMT, t_in); return (n); } t_out = gmt_M_memory (GMT, NULL, n_out + n, double); x_tmp = gmt_M_memory (GMT, NULL, n_out + n, double); y_tmp = gmt_M_memory (GMT, NULL, n_out + n, double); flag = gmt_M_memory (GMT, NULL, n_out + n, bool); /* Create equidistant output points */ ds = t_in[n-1] / (n_out-1); t_next = ds; t_out[0] = 0.0; flag[0] = true; for (i = j = 1; i < n_out; i++) { if (j < n && t_in[j] < t_next) { t_out[i] = t_in[j]; flag[i] = true; j++; n_out++; } else { t_out[i] = t_next; t_next += ds; } } t_out[n_out-1] = t_in[n-1]; if (t_out[n_out-1] == t_out[n_out-2]) n_out--; flag[n_out-1] = true; if (gmt_intpol (GMT, t_in, x, NULL, n, n_out, t_out, x_tmp, 0.0, stype)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT internal error\n"); gmt_M_free (GMT, t_in); gmt_M_free (GMT, t_out); gmt_M_free (GMT, flag); gmt_M_free (GMT, x_tmp); gmt_M_free (GMT, y_tmp); return -GMT_RUNTIME_ERROR; } if (gmt_intpol (GMT, t_in, y, NULL, n, n_out, t_out, y_tmp, 0.0, stype)) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT internal error\n"); gmt_M_free (GMT, t_in); gmt_M_free (GMT, t_out); gmt_M_free (GMT, flag); gmt_M_free (GMT, x_tmp); gmt_M_free (GMT, y_tmp); return -GMT_RUNTIME_ERROR; } /* Make sure interpolated function is bounded on each segment interval */ i = 0; while (i < (n_out - 1)) { j = i + 1; while (j < n_out && !flag[j]) j++; x0 = x_tmp[i]; x1 = x_tmp[j]; if (x0 > x1) gmt_M_double_swap (x0, x1); y0 = y_tmp[i]; y1 = y_tmp[j]; if (y0 > y1) gmt_M_double_swap (y0, y1); for (k = i + 1; k < j; k++) { if (x_tmp[k] < x0) x_tmp[k] = x0 + 1.0e-10; else if (x_tmp[k] > x1) x_tmp[k] = x1 - 1.0e-10; if (y_tmp[k] < y0) y_tmp[k] = y0 + 1.0e-10; else if (y_tmp[k] > y1) y_tmp[k] = y1 - 1.0e-10; } i = j; } /* Replace original coordinates */ gmt_M_free (GMT, x); gmt_M_free (GMT, y); *x_in = x_tmp; *y_in = y_tmp; gmt_M_free (GMT, t_in); gmt_M_free (GMT, t_out); gmt_M_free (GMT, flag); return (n_out); } /*! . */ GMT_LOCAL uint64_t gmtsupport_splice_contour (struct GMT_CTRL *GMT, double **x, double **y, uint64_t n, double *x2, double *y2, uint64_t n2) { /* Basically does a "tail -r" on the array x,y and append arrays x2/y2 */ uint64_t i, j, m; double *x1, *y1; if (n2 < 2) return (n); /* Nothing to be done when second piece < 2 points */ m = n + n2 - 1; /* Total length since one point is shared */ /* Make more space */ x1 = *x; y1 = *y; x1 = gmt_M_memory (GMT, x1, m, double); y1 = gmt_M_memory (GMT, y1, m, double); /* Move first piece to the back */ for (i = m-1, j = n; j > 0; j--, i--) { x1[i] = x1[j-1]; y1[i] = y1[j-1]; } /* Put second piece, in reverse, in the front */ for (i = n2-2, j = 1; j < n2; j++, i--) { x1[i] = x2[j]; y1[i] = y2[j]; } *x = x1; *y = y1; return (m); } /*! . */ GMT_LOCAL void gmtsupport_orient_contour (struct GMT_GRID *G, double *x, double *y, uint64_t n, int orient) { /* Determine handedness of the contour and if opposite of orient reverse the contour */ int side[2], z_dir, k, k2; bool reverse; uint64_t i, j, ij_ul, ij_ur, ij_ll, ij_lr; double fx[2], fy[2], dx, dy; struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (G->header); if (orient == 0) return; /* Nothing to be done when no orientation specified */ if (n < 2) return; /* Cannot work on a single point */ for (k = 0; k < 2; k++) { /* Calculate fractional node numbers from left/top */ fx[k] = (x[k] - G->header->wesn[XLO]) * HH->r_inc[GMT_X] - G->header->xy_off; fy[k] = (G->header->wesn[YHI] - y[k]) * HH->r_inc[GMT_Y] - G->header->xy_off; } /* Get(i,j) of the lower left node in the rectangle containing this contour segment. We use the average x and y coordinate for this to avoid any round-off involved in working on a single coordinate. The average coordinate should always be inside the rectangle and hence the floor/ceil operators will yield the LL node. */ i = lrint (floor (0.5 * (fx[0] + fx[1]))); j = lrint (ceil (0.5 * (fy[0] + fy[1]))); ij_ll = gmt_M_ijp (G->header, j, i); /* lower left corner */ ij_lr = gmt_M_ijp (G->header, j, i+1); /* lower right corner */ ij_ul = gmt_M_ijp (G->header, j-1, i); /* upper left corner */ ij_ur = gmt_M_ijp (G->header, j-1, i+1); /* upper right corner */ for (k = 0; k < 2; k++) { /* Determine which edge the contour points lie on (0-3) */ /* We KNOW that for each k, either x[k] or y[k] lies EXACTLY on a gridline. This is used * to deal with the inevitable round-off that places points slightly off the gridline. We * pick the coordinate closest to the gridline as the one that should be exactly on the gridline */ k2 = 1 - k; /* The other point */ dx = fmod (fx[k], 1.0); if (dx > 0.5) dx = 1.0 - dx; /* Fraction to closest vertical gridline */ dy = fmod (fy[k], 1.0); if (dy > 0.5) dy = 1.0 - dy; /* Fraction to closest horizontal gridline */ if (dx < dy) /* Point is on a vertical grid line (left [3] or right [1]) */ side[k] = (fx[k] < fx[k2]) ? 3 : 1; /* Simply check order of fx to determine which it is */ else /* Point must be on horizontal grid line (top [2] or bottom [0]) */ side[k] = (fy[k] > fy[k2]) ? 0 : 2; /* Same for fy */ } switch (side[0]) { /* Entry side: check heights of corner points.*/ /* if point to the right of the line is higher z_dir = +1 else -1 */ case 0: /* Bottom: check heights of lower left and lower right nodes */ z_dir = (G->data[ij_lr] > G->data[ij_ll]) ? +1 : -1; break; case 1: /* Right */ z_dir = (G->data[ij_ur] > G->data[ij_lr]) ? +1 : -1; break; case 2: /* Top */ z_dir = (G->data[ij_ul] > G->data[ij_ur]) ? +1 : -1; break; default:/* Left */ z_dir = (G->data[ij_ll] > G->data[ij_ul]) ? +1 : -1; break; } reverse = (z_dir != orient); if (reverse) { /* Must reverse order of contour */ for (i = 0, j = n-1; i < n/2; i++, j--) { gmt_M_double_swap (x[i], x[j]); gmt_M_double_swap (y[i], y[j]); } } } /*! . */ GMT_LOCAL bool gmtsupport_label_is_OK (struct GMT_CTRL *GMT, struct GMT_LABEL *L, char *this_label, char *label, double this_dist, double this_value_dist, uint64_t xl, uint64_t fj, struct GMT_CONTOUR *G) { /* Determines if the proposed label passes various tests. Return true if we should go ahead and add this label to the list */ bool label_OK = true; uint64_t seg, k; char format[GMT_LEN256] = {""}; struct GMT_CONTOUR_LINE *S = NULL; if (G->isolate) { /* Must determine if the proposed label is within radius distance of any other label already accepted */ for (seg = 0; seg < G->n_segments; seg++) { /* Previously processed labels */ S = G->segment[seg]; /* Pointer to current segment */ for (k = 0; k < S->n_labels; k++) if (hypot (L->x - S->L[k].x, L->y - S->L[k].y) < G->label_isolation) return (false); } /* Also check labels for current segment */ for (k = 0; k < G->n_label; k++) if (hypot (L->x - G->L[k]->x, L->y - G->L[k]->y) < G->label_isolation) return (false); } switch (G->label_type) { case GMT_LABEL_IS_NONE: if (label && label[0]) strcpy (this_label, label); else label_OK = false; break; case GMT_LABEL_IS_CONSTANT: case GMT_LABEL_IS_HEADER: if (G->label[0]) strcpy (this_label, G->label); else label_OK = false; break; case GMT_LABEL_IS_PDIST: if (G->spacing) { /* Distances are even so use special contour format */ gmt_get_format (GMT, this_dist * GMT->session.u2u[GMT_INCH][G->dist_unit], G->unit, NULL, format); gmt_sprintf_float (GMT, this_label, format, this_dist * GMT->session.u2u[GMT_INCH][G->dist_unit]); } else { gmt_sprintf_float (GMT, this_label, GMT->current.setting.format_float_map, this_dist * GMT->session.u2u[GMT_INCH][G->dist_unit]); } break; case GMT_LABEL_IS_MDIST: gmt_sprintf_float (GMT, this_label, GMT->current.setting.format_float_map, this_value_dist); break; case GMT_LABEL_IS_FFILE: if (G->f_label[fj] && G->f_label[fj][0]) strcpy (this_label, G->f_label[fj]); else label_OK = false; break; case GMT_LABEL_IS_XFILE: if (G->X->table[0]->segment[xl]->label && G->X->table[0]->segment[xl]->label[0]) strcpy (this_label, G->X->table[0]->segment[xl]->label); else label_OK = false; break; case GMT_LABEL_IS_SEG: sprintf (this_label, "%" PRIu64, (GMT->current.io.status & GMT_IO_SEGMENT_HEADER) ? GMT->current.io.seg_no - 1 : GMT->current.io.seg_no); break; case GMT_LABEL_IS_FSEG: sprintf (this_label, "%d/%" PRIu64, GMT->current.io.tbl_no, (GMT->current.io.status & GMT_IO_SEGMENT_HEADER) ? GMT->current.io.seg_no - 1 : GMT->current.io.seg_no); break; default: /* Should not happen... */ GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT internal error - no label set\n"); this_label[0] = '\0'; return false; break; } return (label_OK); } /*! . */ GMT_LOCAL void gmtsupport_place_label (struct GMT_CTRL *GMT, struct GMT_LABEL *L, char *txt, struct GMT_CONTOUR *G, bool use_unit, size_t extra) { /* Allocates needed space and copies in the label */ size_t n, m = 0; if (use_unit && G->unit[0]) m = strlen (G->unit); n = strlen (txt) + 1 + m + extra; if (G->prefix[0]) { /* Must prepend the prefix string */ n += strlen (G->prefix) + 1; L->label = gmt_M_memory (GMT, NULL, n, char); sprintf (L->label, "%s%s", G->prefix, txt); } else { L->label = gmt_M_memory (GMT, NULL, n, char); strcpy (L->label, txt); } if (use_unit && G->unit[0]) { /* Append a unit string */ strcat (L->label, G->unit); } gmt_M_memcpy (L->rgb, G->font_label.fill.rgb, 4, double); /* Remember color in case it varies */ } /*! . */ GMT_LOCAL void gmtsupport_hold_contour_sub (struct GMT_CTRL *GMT, double **xxx, double **yyy, uint64_t nn, double zval, char *label, char ctype, double cangle, bool closed, bool contour, struct GMT_CONTOUR *G) { /* The xxx, yyy are expected to be projected x/y inches */ uint64_t i, j, start = 0; size_t n_alloc = GMT_SMALL_CHUNK; double *track_dist = NULL, *map_dist = NULL, *value_dist = NULL, *radii = NULL, *xx = NULL, *yy = NULL; double dx, dy, width, f, this_dist, step, stept, this_value_dist, lon[2], lat[2]; struct GMT_LABEL *new_label = NULL; char this_label[GMT_BUFSIZ] = {""}; if (nn < 2) return; xx = *xxx; yy = *yyy; G->n_label = 0; /* OK, line is long enough to be added to array of lines */ if (ctype == 'A' || ctype == 'a') { /* Annotated contours, must find label placement */ /* Calculate distance along contour and store in track_dist array */ if (G->dist_kind == 1) gmt_xy_to_geo (GMT, &lon[1], &lat[1], xx[0], yy[0]); map_dist = gmt_M_memory (GMT, NULL, nn, double); /* Distances on map in inches */ track_dist = gmt_M_memory (GMT, NULL, nn, double); /* May be km ,degrees or whatever */ value_dist = gmt_M_memory (GMT, NULL, nn, double); /* May be km ,degrees or whatever */ radii = gmt_M_memory (GMT, NULL, nn, double); /* Radius of curvature, in inches */ /* We will calculate the radii of curvature at all points. By default we don't care and * will place labels at whatever distance we end up with. However, if the user has asked * for a minimum limit on the radius of curvature [Default 0] we do not want to place labels * at those sections where the curvature is large. Since labels are placed according to * distance along track, the way we deal with this is to set distance increments to zero * where curvature is large: that way, there is no increase in distance over those patches * and the machinery for determining when we exceed the next label distance will not kick * in until after curvature drops and increments are again nonzero. This procedure only * applies to the algorithms based on distance along track. */ gmtsupport_get_radii_of_curvature (xx, yy, nn, radii); map_dist[0] = track_dist[0] = value_dist[0] = 0.0; /* Unnecessary, just so we understand the logic */ for (i = 1; i < nn; i++) { /* Distance from xy */ dx = xx[i] - xx[i-1]; if (gmt_M_is_geographic (GMT, GMT_IN) && GMT->current.map.is_world && fabs (dx) > (width = gmt_half_map_width (GMT, yy[i-1]))) { width *= 2.0; dx = copysign (width - fabs (dx), -dx); if (xx[i] < width) xx[i-1] -= width; else xx[i-1] += width; } dy = yy[i] - yy[i-1]; step = stept = hypot (dx, dy); map_dist[i] = map_dist[i-1] + step; if (G->dist_kind == 1 || G->label_type == GMT_LABEL_IS_MDIST) { lon[0] = lon[1]; lat[0] = lat[1]; gmt_xy_to_geo (GMT, &lon[1], &lat[1], xx[i], yy[i]); if (G->dist_kind == 1) step = gmtlib_distance_type (GMT, lon[0], lat[0], lon[1], lat[1], GMT_CONT_DIST); if (G->label_type == GMT_LABEL_IS_MDIST) stept = gmtlib_distance_type (GMT, lon[0], lat[0], lon[1], lat[1], GMT_LABEL_DIST); } if (radii[i] < G->min_radius) step = stept = 0.0; /* If curvature is too great we simply don't add up distances */ track_dist[i] = track_dist[i-1] + step; value_dist[i] = value_dist[i-1] + stept; } gmt_M_free (GMT, radii); /* G->L array is only used so we can later sort labels based on distance along track. Once * GMT_contlabel_draw has been called we will free up the memory as the labels are kept in * the linked list starting at G->anchor. */ G->L = gmt_M_memory (GMT, NULL, n_alloc, struct GMT_LABEL *); if (G->spacing) { /* Place labels based on distance along contours */ double last_label_dist, dist_offset, dist; dist_offset = (closed && G->dist_kind == 0) ? (1.0 - G->label_dist_frac) * G->label_dist_spacing : 0.0; /* Label closed contours longer than frac of dist_spacing */ last_label_dist = 0.0; i = 1; while (i < nn) { dist = track_dist[i] + dist_offset - last_label_dist; if (dist > G->label_dist_spacing) { /* Time for label */ new_label = gmt_M_memory (GMT, NULL, 1, struct GMT_LABEL); f = (dist - G->label_dist_spacing) / (track_dist[i] - track_dist[i-1]); if (f < 0.5) { new_label->x = xx[i] - f * (xx[i] - xx[i-1]); new_label->y = yy[i] - f * (yy[i] - yy[i-1]); new_label->dist = map_dist[i] - f * (map_dist[i] - map_dist[i-1]); this_value_dist = value_dist[i] - f * (value_dist[i] - value_dist[i-1]); } else { f = 1.0 - f; new_label->x = xx[i-1] + f * (xx[i] - xx[i-1]); new_label->y = yy[i-1] + f * (yy[i] - yy[i-1]); new_label->dist = map_dist[i-1] + f * (map_dist[i] - map_dist[i-1]); this_value_dist = value_dist[i-1] + f * (value_dist[i] - value_dist[i-1]); } this_dist = G->label_dist_spacing - dist_offset + last_label_dist; if (gmtsupport_label_is_OK (GMT, new_label, this_label, label, this_dist, this_value_dist, 0, 0, G)) { gmtsupport_place_label (GMT, new_label, this_label, G, !(G->label_type == GMT_LABEL_IS_NONE || G->label_type == GMT_LABEL_IS_PDIST), 0); new_label->node = i - 1; gmtsupport_contlabel_angle (GMT, xx, yy, i - 1, i, cangle, nn, contour, new_label, G); G->L[G->n_label++] = new_label; if (G->n_label == n_alloc) { size_t old_n_alloc = n_alloc; n_alloc <<= 1; G->L = gmt_M_memory (GMT, G->L, n_alloc, struct GMT_LABEL *); gmt_M_memset (&(G->L[old_n_alloc]), n_alloc - old_n_alloc, struct GMT_LABEL *); /* Set to NULL */ } } else /* All in vain... */ gmt_M_free (GMT, new_label); dist_offset = 0.0; last_label_dist = this_dist; } else /* Go to next point in line */ i++; } if (G->n_label == 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your -Gd|D option produced no contour labels for z = %g\n", zval); } if (G->number) { /* Place prescribed number of labels evenly along contours */ uint64_t nc; int e_val = 0; double dist, last_dist; last_dist = (G->n_cont > 1) ? -map_dist[nn-1] / (G->n_cont - 1) : -0.5 * map_dist[nn-1]; nc = (map_dist[nn-1] > G->min_dist) ? G->n_cont : 0; for (i = j = 0; i < nc; i++) { new_label = gmt_M_memory (GMT, NULL, 1, struct GMT_LABEL); if (G->number_placement && !closed) { e_val = G->number_placement; if (G->number_placement == -1 && G->n_cont == 1) { /* Label justified with start of segment */ f = d_atan2d (xx[0] - xx[1], yy[0] - yy[1]) + 180.0; /* 0-360 */ G->end_just[0] = (f >= 90.0 && f <= 270) ? 7 : 5; } else if (G->number_placement == +1 && G->n_cont == 1) { /* Label justified with end of segment */ f = d_atan2d (xx[nn-1] - xx[nn-2], yy[nn-1] - yy[nn-2]) + 180.0; /* 0-360 */ G->end_just[1] = (f >= 90.0 && f <= 270) ? 7 : 5; } else if (G->number_placement && G->n_cont > 1) /* One of the end labels */ e_val = (i == 0) ? -1 : +1; dist = (G->n_cont > 1) ? i * track_dist[nn-1] / (G->n_cont - 1) : 0.5 * (G->number_placement + 1.0) * track_dist[nn-1]; this_value_dist = (G->n_cont > 1) ? i * value_dist[nn-1] / (G->n_cont - 1) : 0.5 * (G->number_placement + 1.0) * value_dist[nn-1]; } else { dist = (i + 1 - 0.5 * closed) * track_dist[nn-1] / (G->n_cont + 1 - closed); this_value_dist = (i + 1 - 0.5 * closed) * value_dist[nn-1] / (G->n_cont + 1 - closed); } while (j < nn && track_dist[j] < dist) j++; if (j == nn) j--; /* Safety precaution */ f = ((j == 0) ? 1.0 : (dist - track_dist[j-1]) / (track_dist[j] - track_dist[j-1])); if (f < 0.5) { /* Pick the smallest fraction to minimize Taylor shortcomings */ new_label->x = xx[j-1] + f * (xx[j] - xx[j-1]); new_label->y = yy[j-1] + f * (yy[j] - yy[j-1]); new_label->dist = map_dist[j-1] + f * (map_dist[j] - map_dist[j-1]); } else { f = 1.0 - f; new_label->x = (j == 0) ? xx[0] : xx[j] - f * (xx[j] - xx[j-1]); new_label->y = (j == 0) ? yy[0] : yy[j] - f * (yy[j] - yy[j-1]); new_label->dist = (j == 0) ? 0.0 : map_dist[j] - f * (map_dist[j] - map_dist[j-1]); } if ((new_label->dist - last_dist) >= G->min_dist) { /* OK to accept this label */ this_dist = dist; if (gmtsupport_label_is_OK (GMT, new_label, this_label, label, this_dist, this_value_dist, 0, 0, G)) { size_t extra = (G->crossect) ? strlen (G->crossect_tag[i]) + 1 : 0; /* Need to increase allocated space */ gmtsupport_place_label (GMT, new_label, this_label, G, !(G->label_type == GMT_LABEL_IS_NONE), extra); if (G->crossect) { /* Special crossection mode */ if (!strcmp (new_label->label, "N/A")) /* Override the N/A lack of label identifier */ strcpy (new_label->label, G->crossect_tag[i]); else /* Append tag to the label */ strcat (new_label->label, G->crossect_tag[i]); } new_label->node = (j == 0) ? 0 : j - 1; gmtsupport_contlabel_angle (GMT, xx, yy, new_label->node, j, cangle, nn, contour, new_label, G); if (G->number_placement) new_label->end = e_val; G->L[G->n_label++] = new_label; if (G->n_label == n_alloc) { size_t old_n_alloc = n_alloc; n_alloc <<= 1; G->L = gmt_M_memory (GMT, G->L, n_alloc, struct GMT_LABEL *); gmt_M_memset (&(G->L[old_n_alloc]), n_alloc - old_n_alloc, struct GMT_LABEL *); /* Set to NULL */ } last_dist = new_label->dist; } else /* Label had no text */ gmt_M_free (GMT, new_label); } else /* All in vain... */ gmt_M_free (GMT, new_label); } if (G->n_label == 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your -Gn|N option produced no contour labels for z = %g\n", zval); } if (G->crossing) { /* Determine label positions based on crossing lines */ uint64_t left, right, line_no; struct GMT_DATASEGMENT *S = NULL; gmt_init_track (GMT, yy, nn, &(G->ylist)); for (line_no = 0; line_no < G->X->n_segments; line_no++) { /* For each of the crossing lines */ S = G->X->table[0]->segment[line_no]; /* Current segment */ gmt_init_track (GMT, S->data[GMT_Y], S->n_rows, &(G->ylist_XP)); G->nx = (unsigned int)gmt_crossover (GMT, S->data[GMT_X], S->data[GMT_Y], NULL, G->ylist_XP, S->n_rows, xx, yy, NULL, G->ylist, nn, false, gmt_M_is_geographic (GMT, GMT_IN), &G->XC); gmt_M_free (GMT, G->ylist_XP); if (G->nx == 0) continue; /* OK, we found intersections for labels */ for (i = 0; i < (uint64_t)G->nx; i++) { left = lrint (floor (G->XC.xnode[1][i])); right = lrint (ceil (G->XC.xnode[1][i])); new_label = gmt_M_memory (GMT, NULL, 1, struct GMT_LABEL); new_label->x = G->XC.x[i]; new_label->y = G->XC.y[i]; new_label->node = left; new_label->dist = track_dist[left]; f = G->XC.xnode[1][i] - left; if (f < 0.5) { this_dist = track_dist[left] + f * (track_dist[right] - track_dist[left]); new_label->dist = map_dist[left] + f * (map_dist[right] - map_dist[left]); this_value_dist = value_dist[left] + f * (value_dist[right] - value_dist[left]); } else { f = 1.0 - f; this_dist = track_dist[right] - f * (track_dist[right] - track_dist[left]); new_label->dist = map_dist[right] - f * (map_dist[right] - map_dist[left]); this_value_dist = value_dist[right] - f * (value_dist[right] - value_dist[left]); } if (gmtsupport_label_is_OK (GMT, new_label, this_label, label, this_dist, this_value_dist, line_no, 0, G)) { gmtsupport_place_label (GMT, new_label, this_label, G, !(G->label_type == GMT_LABEL_IS_NONE), 0); gmtsupport_contlabel_angle (GMT, xx, yy, left, right, cangle, nn, contour, new_label, G); G->L[G->n_label++] = new_label; if (G->n_label == n_alloc) { size_t old_n_alloc = n_alloc; n_alloc <<= 1; G->L = gmt_M_memory (GMT, G->L, n_alloc, struct GMT_LABEL *); gmt_M_memset (&(G->L[old_n_alloc]), n_alloc - old_n_alloc, struct GMT_LABEL *); /* Set to NULL */ } } else /* All in vain... */ gmt_M_free (GMT, new_label); } gmt_x_free (GMT, &G->XC); } gmt_M_free (GMT, G->ylist); if (G->n_label == 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your -Gx|X|l|L option produced no contour labels for z = %g\n", zval); } if (G->fixed) { /* Prescribed point locations for labels that match points in input records */ double dist, min_dist; for (j = 0; j < (uint64_t)G->f_n; j++) { /* Loop over fixed point list */ min_dist = DBL_MAX; for (i = 0; i < nn; i++) { /* Loop over input line/contour */ if ((dist = hypot (xx[i] - G->f_xy[0][j], yy[i] - G->f_xy[1][j])) < min_dist) { /* Better fit */ min_dist = dist; start = i; } } if (min_dist < G->slop) { /* Closest point within tolerance */ new_label = gmt_M_memory (GMT, NULL, 1, struct GMT_LABEL); new_label->x = xx[start]; new_label->y = yy[start]; new_label->node = start; new_label->dist = track_dist[start]; this_dist = track_dist[start]; new_label->dist = map_dist[start]; this_value_dist = value_dist[start]; if (gmtsupport_label_is_OK (GMT, new_label, this_label, label, this_dist, this_value_dist, 0, j, G)) { gmtsupport_place_label (GMT, new_label, this_label, G, !(G->label_type == GMT_LABEL_IS_NONE), 0); gmtsupport_contlabel_angle (GMT, xx, yy, start, start, cangle, nn, contour, new_label, G); G->L[G->n_label++] = new_label; if (G->n_label == n_alloc) { size_t old_n_alloc = n_alloc; n_alloc <<= 1; G->L = gmt_M_memory (GMT, G->L, n_alloc, struct GMT_LABEL *); gmt_M_memset (&(G->L[old_n_alloc]), n_alloc - old_n_alloc, struct GMT_LABEL *); /* Set to NULL */ } } else /* All in vain... */ gmt_M_free (GMT, new_label); } } if (G->n_label == 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your -Gf option produced no contour labels for z = %g\n", zval); } gmtsupport_contlabel_fixpath (GMT, &xx, &yy, map_dist, &nn, G); /* Inserts the label x,y into path */ gmtsupport_contlabel_addpath (GMT, xx, yy, nn, zval, label, true, G); /* Appends this path and the labels to list */ gmt_M_free (GMT, track_dist); gmt_M_free (GMT, map_dist); gmt_M_free (GMT, value_dist); /* Free label structure since info is now copied to segment labels */ for (i = 0; i < (uint64_t)G->n_label; i++) { gmt_M_free (GMT, G->L[i]->label); gmt_M_free (GMT, G->L[i]); } gmt_M_free (GMT, G->L); } else { /* just one line, no holes for labels */ gmtsupport_contlabel_addpath (GMT, xx, yy, nn, zval, label, false, G); /* Appends this path to list */ } *xxx = xx; *yyy = yy; } /*! . */ GMT_LOCAL void gmtsupport_add_decoration (struct GMT_CTRL *GMT, struct GMT_DATASEGMENT *S, struct GMT_LABEL *L, struct GMT_DECORATE *G) { /* Add a symbol location to the growing segment */ struct GMT_DATASEGMENT_HIDDEN *SH = gmt_get_DS_hidden (S); if (S->n_rows == SH->n_alloc) { /* Need more memory for the segment */ uint64_t col; SH->n_alloc += GMT_SMALL_CHUNK; for (col = 0; col < S->n_columns; col++) S->data[col] = gmt_M_memory (GMT, S->data[col], SH->n_alloc, double); S->text = gmt_M_memory (GMT, S->text, SH->n_alloc, char *); SH->alloc_mode = GMT_ALLOC_INTERNALLY; } /* Deal with any justifications or nudging */ if (G->nudge_flag) { /* Must adjust point a bit */ double s = 0.0, c = 1.0, sign = 1.0; if (G->nudge_flag == 2) sincosd (L->line_angle, &s, &c); /* If N+1 or N-1 is used we want positive x nudge to extend away from end point */ sign = (G->number_placement) ? (double)L->end : 1.0; L->x += sign * (G->nudge[GMT_X] * c - G->nudge[GMT_Y] * s); L->y += sign * (G->nudge[GMT_X] * s + G->nudge[GMT_Y] * c); } /* Build record with Cartesian (hence col = GMT_Z) "x y angle symbol" since we are using -Jx1i */ S->data[GMT_X][S->n_rows] = L->x; S->data[GMT_Y][S->n_rows] = L->y; S->data[GMT_Z][S->n_rows] = gmt_M_to_inch (GMT, G->size); S->data[3][S->n_rows] = L->angle; /* Change this in inches internally instead of string */ S->text[S->n_rows++] = strdup (G->symbol_code); } /*! . */ GMT_LOCAL void gmtsupport_decorated_line_sub (struct GMT_CTRL *GMT, double *xx, double *yy, uint64_t nn, struct GMT_DECORATE *G, struct GMT_DATASET *D, uint64_t seg) { /* The xxx, yyy are expected to be projected x/y inches. * This function is modelled after gmtsupport_hold_contour_sub but tweaked to deal with * the placement of psxy-symbols rather that text labels. This is in most regards * simpler than placing text so many lines related to text have been yanked. There * is the assumption that the symbols will be filled so we make no attempt to clip * the lines as we do for contours/quated lines. It would be very complicated to * determine that clip path if the symbol is a star, for instance. So if fill is * not given then we will see the line beneath the symbol. */ uint64_t i, j, start = 0; bool closed; double *track_dist = NULL, *map_dist = NULL, *value_dist = NULL; double dx, dy, width, f, this_dist, step, stept, lon[2], lat[2]; struct GMT_DATASEGMENT *S = D->table[0]->segment[seg]; /* A single segment */ struct GMT_LABEL L; /* Needed to pick up angles */ if (nn < 2) return; /* You, sir, are not a line! */ closed = (nn > 2 && !gmt_polygon_is_open (GMT, xx, yy, nn)); /* true if this is a polygon */ /* Calculate distance along line and store in track_dist array */ if (G->dist_kind == 1) gmt_xy_to_geo (GMT, &lon[1], &lat[1], xx[0], yy[0]); map_dist = gmt_M_memory (GMT, NULL, nn, double); /* Distances on map in inches */ track_dist = gmt_M_memory (GMT, NULL, nn, double); /* May be km, degrees or whatever */ value_dist = gmt_M_memory (GMT, NULL, nn, double); /* May be km, degrees or whatever */ map_dist[0] = track_dist[0] = value_dist[0] = 0.0; /* Unnecessary, just so we understand the logic */ for (i = 1; i < nn; i++) { /* Distance from xy in plot distances (inch) */ dx = xx[i] - xx[i-1]; if (gmt_M_is_geographic (GMT, GMT_IN) && GMT->current.map.is_world && fabs (dx) > (width = gmt_half_map_width (GMT, yy[i-1]))) { width *= 2.0; dx = copysign (width - fabs (dx), -dx); if (xx[i] < width) xx[i-1] -= width; else xx[i-1] += width; } dy = yy[i] - yy[i-1]; step = stept = hypot (dx, dy); /* Initially these steps are in inches */ map_dist[i] = map_dist[i-1] + step; if (G->dist_kind == 1) { /* Wanted spacing in map distance units */ lon[0] = lon[1]; lat[0] = lat[1]; gmt_xy_to_geo (GMT, &lon[1], &lat[1], xx[i], yy[i]); if (G->dist_kind == 1) step = gmtlib_distance_type (GMT, lon[0], lat[0], lon[1], lat[1], GMT_CONT_DIST); } track_dist[i] = track_dist[i-1] + step; value_dist[i] = value_dist[i-1] + stept; } if (G->spacing) { /* Place symbols based on distance along lines */ double last_label_dist, dist_offset, dist; dist_offset = (closed && G->dist_kind == 0) ? (1.0 - G->symbol_dist_frac) * G->symbol_dist_spacing : 0.0; /* Only place symbols on closed contours longer than frac of dist_spacing */ last_label_dist = 0.0; i = 1; while (i < nn) { dist = track_dist[i] + dist_offset - last_label_dist; if (dist > G->symbol_dist_spacing) { /* Time for symbol */ f = (dist - G->symbol_dist_spacing) / (track_dist[i] - track_dist[i-1]); if (f < 0.5) { L.x = xx[i] - f * (xx[i] - xx[i-1]); L.y = yy[i] - f * (yy[i] - yy[i-1]); } else { f = 1.0 - f; L.x = xx[i-1] + f * (xx[i] - xx[i-1]); L.y = yy[i-1] + f * (yy[i] - yy[i-1]); } this_dist = G->symbol_dist_spacing - dist_offset + last_label_dist; gmtsupport_decorated_angle (GMT, xx, yy, i - 1, i, G->symbol_angle, nn, false, &L, G); gmtsupport_add_decoration (GMT, S, &L, G); dist_offset = 0.0; last_label_dist = this_dist; } else /* Go to next point in line */ i++; } } else if (G->number) { /* Place prescribed number of symbols evenly along lines */ uint64_t nc; double dist; nc = (map_dist[nn-1] > G->min_dist) ? G->n_cont : 0; L.end = 0; for (i = j = 0; i < nc; i++) { if (G->number_placement && !closed) { dist = (G->n_cont > 1) ? i * track_dist[nn-1] / (G->n_cont - 1) : 0.5 * (G->number_placement + 1.0) * track_dist[nn-1]; L.end = (G->number_placement && G->n_cont > 1) ? ((i == 0) ? -1 : +1) : G->number_placement; } else dist = (i + 1 - 0.5 * closed) * track_dist[nn-1] / (G->n_cont + 1 - closed); while (j < nn && track_dist[j] < dist) j++; if (j == nn) j--; /* Safety precaution */ f = ((j == 0) ? 1.0 : (dist - track_dist[j-1]) / (track_dist[j] - track_dist[j-1])); if (f < 0.5) { /* Pick the smallest fraction to minimize Taylor shortcomings */ L.x = xx[j-1] + f * (xx[j] - xx[j-1]); L.y = yy[j-1] + f * (yy[j] - yy[j-1]); } else { f = 1.0 - f; L.x = (j == 0) ? xx[0] : xx[j] - f * (xx[j] - xx[j-1]); L.y = (j == 0) ? yy[0] : yy[j] - f * (yy[j] - yy[j-1]); } L.node = (j == 0) ? 0 : j - 1; gmtsupport_decorated_angle (GMT, xx, yy, L.node, j, G->symbol_angle, nn, false, &L, G); gmtsupport_add_decoration (GMT, S, &L, G); } } else if (G->crossing) { /* Determine label positions based on crossing lines */ uint64_t left, right, line_no; struct GMT_DATASEGMENT *Sd = NULL; gmt_init_track (GMT, yy, nn, &(G->ylist)); for (line_no = 0; line_no < G->X->n_segments; line_no++) { /* For each of the crossing lines */ Sd = G->X->table[0]->segment[line_no]; /* Current segment */ gmt_init_track (GMT, Sd->data[GMT_Y], Sd->n_rows, &(G->ylist_XP)); G->nx = (unsigned int)gmt_crossover (GMT, Sd->data[GMT_X], Sd->data[GMT_Y], NULL, G->ylist_XP, Sd->n_rows, xx, yy, NULL, G->ylist, nn, false, gmt_M_is_geographic (GMT, GMT_IN), &G->XC); gmt_M_free (GMT, G->ylist_XP); if (G->nx == 0) continue; /* OK, we found intersections for labels */ for (i = 0; i < (uint64_t)G->nx; i++) { left = lrint (floor (G->XC.xnode[1][i])); right = lrint (ceil (G->XC.xnode[1][i])); L.x = G->XC.x[i]; L.y = G->XC.y[i]; gmtsupport_decorated_angle (GMT, xx, yy, left, right, G->symbol_angle, nn, false, &L, G); gmtsupport_add_decoration (GMT, S, &L, G); } gmt_x_free (GMT, &G->XC); } gmt_M_free (GMT, G->ylist); } else if (G->fixed) { /* Prescribed point locations for labels that match points in input records */ double dist, min_dist; for (j = 0; j < (uint64_t)G->f_n; j++) { /* Loop over fixed point list */ min_dist = DBL_MAX; for (i = 0; i < nn; i++) { /* Loop over input line/contour */ if ((dist = hypot (xx[i] - G->f_xy[0][j], yy[i] - G->f_xy[1][j])) < min_dist) { /* Better fit */ min_dist = dist; start = i; } } if (min_dist < G->slop) { /* Closest point within tolerance */ L.x = xx[start]; L.y = yy[start]; gmtsupport_decorated_angle (GMT, xx, yy, start, start, G->symbol_angle, nn, false, &L, G); gmtsupport_add_decoration (GMT, S, &L, G); } } } gmt_M_free (GMT, track_dist); gmt_M_free (GMT, map_dist); gmt_M_free (GMT, value_dist); } /*! . */ GMT_LOCAL uint64_t gmtsupport_getprevpoint (double plon, double lon[], uint64_t n, uint64_t this_p) { /* Return the previous point that does NOT equal plon */ uint64_t ip = (this_p == 0) ? n - 2 : this_p - 1; /* Previous point (-2 because last is a duplicate of first) */ while (doubleAlmostEqualZero (plon, lon[ip]) || doubleAlmostEqual (fabs(plon - lon[ip]), 360.0)) { /* Same as plon */ if (ip == 0) ip = n - 2; else ip--; } return (ip); } /*! . */ static inline bool gmtsupport_same_longitude (double a, double b) { /* return true if a and b are the same longitude */ while (a < 0.0) a += 360.0; while (a > 360.0) a -= 360.0; while (b < 0.0) b += 360.0; while (b > 360.0) b -= 360.0; return doubleAlmostEqualZero (a, b); } #define GMT_SAME_LATITUDE(A,B) (doubleAlmostEqualZero (A,B)) /* A and B are the same latitude */ /*! . */ GMT_LOCAL int gmtsupport_inonout_sphpol_count (double plon, double plat, const struct GMT_DATASEGMENT *P, unsigned int count[]) { /* Case of a polar cap */ uint64_t i, in, ip, prev; int cut; double W, E, S, N, lon, lon1, lon2, dlon, x_lat, dx1, dx2; /* Draw meridian through P and count all the crossings with the line segments making up the polar cap S */ gmt_M_memset (count, 2, unsigned int); /* Initialize counts to zero */ for (i = 0; i < P->n_rows - 1; i++) { /* -1, since we know last point repeats the first */ /* Here lon1 and lon2 are the end points (in longitude) of the current line segment in S. There are * four cases to worry about: * 1) lon equals lon1 (i.e., the meridian through lon goes right through lon1) * 2) lon equals lon2 (i.e., the meridian through lon goes right through lon2) * 3) lon lies between lon1 and lon2 and crosses the segment * 4) none of the above * Since we want to obtain either ONE or ZERO intersections per segment we will skip to next * point if case (2) occurs: this avoids counting a crossing twice for consecutive segments. */ if (gmtsupport_same_longitude (plon, P->data[GMT_X][i]) && GMT_SAME_LATITUDE (plat, P->data[GMT_Y][i])) return (GMT_ONEDGE); /* Point is on the perimeter */ in = i + 1; /* Next point index */ /* Next deal with case when the longitude of P goes ~right through the second of the line nodes */ if (gmtsupport_same_longitude (plon, P->data[GMT_X][in])) continue; /* Line goes through the 2nd node - ignore */ lon1 = P->data[GMT_X][i]; /* Copy the first of two longitudes since we may need to mess with them */ lon2 = P->data[GMT_X][in]; /* Copy the second of two longitudes since we may need to mess with them */ if (gmtsupport_same_longitude (plon, lon1)) { /* Line goes through the 1st node */ /* Must check that the two neighboring points are on either side; otherwise it is just a tangent line */ ip = gmtsupport_getprevpoint (plon, P->data[GMT_X], P->n_rows, i); /* Index of previous point != plon */ gmt_M_set_delta_lon (lon1, P->data[GMT_X][ip], dx1); /* Allow for jumps across discontinuous 0 or 180 boundary */ if (dx1 == 0.0) continue; /* Points ip and i forms a meridian, we a tangent line */ gmt_M_set_delta_lon (lon1, lon2, dx2); /* Allow for jumps across discontinuous 0 or 180 boundary */ if (dx1*dx2 > 0.0) continue; /* Both on same side since signs are the same */ cut = (P->data[GMT_Y][i] > plat) ? 0 : 1; /* node is north (0) or south (1) of P */ count[cut]++; prev = ip + 1; /* Always exists because ip is <= n-2 */ /* If prev < i then we have a vertical segment of 2 or more points; prev points to the other end of the segment. * We must then check if our points plat is within that range, meaning the point lies on the segment */ if (prev < i && ((plat <= P->data[GMT_Y][prev] && plat >= P->data[GMT_Y][i]) || (plat <= P->data[GMT_Y][i] && plat >= P->data[GMT_Y][prev]))) return (GMT_ONEDGE); /* P is on segment boundary; we are done */ continue; } /* OK, not exactly on a node, deal with crossing a line */ dlon = lon2 - lon1; if (dlon > 180.0) /* Jumped across Greenwich going westward */ lon2 -= 360.0; else if (dlon < -180.0) /* Jumped across Greenwich going eastward */ lon1 -= 360.0; if (lon1 <= lon2) { /* Segment goes W to E (or N-S) */ W = lon1; E = lon2; } else { /* Segment goes E to W */ W = lon2; E = lon1; } lon = plon; /* Local copy of plon, below adjusted given the segment lon range */ while (lon > W) lon -= 360.0; /* Make sure we rewind way west for starters */ while (lon < W) lon += 360.0; /* Then make sure we wind to inside the lon range or way east */ if (lon > E) continue; /* Not crossing this segment */ if (dlon == 0.0) { /* Special case of N-S segment: does P lie on it? */ if (P->data[GMT_Y][in] < P->data[GMT_Y][i]) { /* Get N and S limits for segment */ S = P->data[GMT_Y][in]; N = P->data[GMT_Y][i]; } else { N = P->data[GMT_Y][in]; S = P->data[GMT_Y][i]; } if (plat < S || plat > N) continue; /* P is not on this segment */ return (1); /* P is on segment boundary; we are done*/ } /* Calculate latitude at intersection */ if (GMT_SAME_LATITUDE (P->data[GMT_Y][i], P->data[GMT_Y][in])) { /* Special cases */ if (GMT_SAME_LATITUDE (plat, P->data[GMT_Y][in])) return (GMT_ONEDGE); /* P is on S boundary */ x_lat = P->data[GMT_Y][i]; } else x_lat = P->data[GMT_Y][i] + ((P->data[GMT_Y][in] - P->data[GMT_Y][i]) / (lon2 - lon1)) * (lon - lon1); if (doubleAlmostEqualZero (x_lat, plat)) return (1); /* P is on S boundary */ cut = (x_lat > plat) ? 0 : 1; /* Cut is north (0) or south (1) of P */ count[cut]++; } return (0); /* This means no special cases were detected that warranted an immediate return */ } /*! . */ GMT_LOCAL unsigned int gmtsupport_inonout_sphpol (struct GMT_CTRL *GMT, double plon, double plat, struct GMT_DATASEGMENT *S) { /* This function is used to see if some point P = (plon, plat) is located inside, outside, or on the boundary of the * spherical polygon S read by GMT_import_table. Note GMT->current.io.skip_duplicates must be true when the polygon * was read so there are NO duplicate (repeated) points. * Returns the following values: * 0: P is outside of S * 1: P is inside of S * 2: P is on boundary of S */ /* Algorithm: * Case 1: The polygon S contains a geographical pole * a) if P is beyond the far latitude then P is outside * b) Draw meridian through P and count intersections: * odd: P is outside; even: P is inside * Case 2: S does not contain a pole * a) If P is outside range of latitudes then P is outside * c) Draw meridian through P and count intersections: * odd: P is inside; even: P is outside * In all cases, we check if P is on the outline of S * Limitation: We assume points are closely spaced so that we can do linear * approximation between successive points in the polygon. */ unsigned int count[2]; struct GMT_DATASEGMENT_HIDDEN *SH = gmt_get_DS_hidden (S); gmt_M_unused(GMT); if (SH->pole) { /* Case 1 of an enclosed polar cap */ if (SH->pole == +1) { /* N polar cap */ if (plat < S->min[GMT_Y]) return (GMT_OUTSIDE); /* South of a N polar cap */ if (plat > SH->lat_limit) return (GMT_INSIDE); /* Clearly inside of a N polar cap */ } else if (SH->pole == -1) { /* S polar cap */ if (plat > S->max[GMT_Y]) return (GMT_OUTSIDE); /* North of a S polar cap */ if (plat < SH->lat_limit) return (GMT_INSIDE); /* Clearly inside of a S polar cap */ } /* Tally up number of intersections between polygon and meridian through P */ if (gmtsupport_inonout_sphpol_count (plon, plat, S, count)) return (GMT_ONEDGE); /* Found P is on S */ if (SH->pole == +1 && count[0] % 2 == 0) return (GMT_INSIDE); if (SH->pole == -1 && count[1] % 2 == 0) return (GMT_INSIDE); return (GMT_OUTSIDE); } /* Here is Case 2. First check latitude range */ if (plat < S->min[GMT_Y] || plat > S->max[GMT_Y]) return (GMT_OUTSIDE); /* Longitudes are tricker and are tested with the tallying of intersections */ if (gmtsupport_inonout_sphpol_count (plon, plat, S, count)) return (GMT_ONEDGE); /* Found P is on S */ if (count[0] % 2) return (GMT_INSIDE); return (GMT_OUTSIDE); /* Nothing triggered the tests; we are outside */ } /*! . */ GMT_LOCAL unsigned int gmtsupport_inonout_sub (struct GMT_CTRL *GMT, double x, double y, struct GMT_DATASEGMENT *S) { /* Front end for both spherical and Cartesian in-on-out functions */ unsigned int side; if (GMT->current.proj.sph_inside) { /* Assumes these are input polygons */ struct GMT_DATASEGMENT_HIDDEN *SH = gmt_get_DS_hidden (S); if (SH->pole) /* 360-degree polar cap, must check fully */ side = gmtsupport_inonout_sphpol (GMT, x, y, S); else { /* See if we are outside range of longitudes for polygon */ while (x > S->min[GMT_X]) x -= 360.0; /* Wind clear of west */ while (x < S->min[GMT_X]) x += 360.0; /* Wind east until inside or beyond east */ if (x > S->max[GMT_X]) return (GMT_OUTSIDE); /* Point outside, no need to assign value */ side = gmtsupport_inonout_sphpol (GMT, x, y, S); } } else { /* Flat Earth case */ if (y < S->min[GMT_Y] || y > S->max[GMT_Y]) return (GMT_OUTSIDE); /* Point outside, no need to assign value */ if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Deal with longitude periodicity */ if (x < S->min[GMT_X]) { x += 360.0; if (x > S->max[GMT_X]) return (GMT_OUTSIDE); /* Point outside, no need to assign value */ } else if (x > S->max[GMT_X]) { x -= 360.0; if (x < S->min[GMT_X]) return (GMT_OUTSIDE); /* Point outside, no need to assign value */ } } else if (x < S->min[GMT_X] || x > S->max[GMT_X]) return (GMT_OUTSIDE); /* Point outside, no need to assign value */ /* It is essential that the longitudes in S are within the min/max limits since this is a Cartesian algorithm */ side = gmt_non_zero_winding (GMT, x, y, S->data[GMT_X], S->data[GMT_Y], S->n_rows); } return (side); } int gmt_signum (double x) { /* The standard sign function */ if (x < 0.0) return -1; if (x > 0.0) return +1; return 0; } #ifdef TRIANGLE_D /* * New gmt_delaunay interface routine that calls the triangulate function * developed by Jonathan Richard Shewchuk, University of California at Berkeley. * Suggested by alert GMT user Alain Coat. You need to get triangle.c and * triangle.h from www.cs.cmu.edu/~quake/triangle.html */ #define REAL double #include "triangle.h" /* Leave link as int**, not uint64_t** */ /*! . */ GMT_LOCAL uint64_t gmtsupport_delaunay_shewchuk (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n, int **link) { /* GMT interface to the triangle package; see above for references. * All that is done is reformatting of parameters and calling the * main triangulate routine. Thanx to Alain Coat for the tip. */ uint64_t i, j; struct triangulateio In, Out, vorOut; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Delaunay triangulation calculated by Jonathan Shewchuk's Triangle [http://www.cs.cmu.edu/~quake/triangle.html]\n"); /* Set everything to 0 and NULL */ gmt_M_memset (&In, 1, struct triangulateio); gmt_M_memset (&Out, 1, struct triangulateio); gmt_M_memset (&vorOut, 1, struct triangulateio); /* Allocate memory for input points */ In.numberofpoints = (int)n; In.pointlist = gmt_M_memory (GMT, NULL, 2 * n, double); /* Copy x,y points to In structure array */ for (i = j = 0; i < n; i++) { In.pointlist[j++] = x_in[i]; In.pointlist[j++] = y_in[i]; } /* Call Jonathan Shewchuk's triangulate algorithm. This is 64-bit safe since * all the structures use 4-byte ints (longs are used internally). The options are * z : Start numbering at zero instead of 1. * I : No iteration numbers * Q : Quiet - suppress all explanations * B : Suppresses output of boundary information * D : Conforming Delaunay: all triangles are truly Delaunay. */ triangulate ("zIQBD", &In, &Out, &vorOut); *link = Out.trianglelist; /* List of node numbers to return via link [NOT ALLOCATED BY gmt_M_memory] */ gmt_M_str_free (Out.pointlist); gmt_M_free (GMT, In.pointlist); return (Out.numberoftriangles); } /*! . */ GMT_LOCAL struct GMT_DATASET * gmtsupport_voronoi_shewchuk (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n64, double *wesn, unsigned int mode) { /* GMT interface to the triangle package; see above for references. * All that is done is reformatting of parameters and calling the * main triangulate routine. Here we return Voronoi information * and package the coordinates of the edges in the output dataset. * The wesn[] array contains the min/max x (or lon) and y (or lat) coordinates. * Since the internal variables in triangle are ints (e.g., edgelist) we use * ints here as well */ /* Currently we only write the edges of a Voronoi cell but we want polygons later. * Info from triangle re Voronoi polygons: "Triangle does not write a list of * the edges adjoining each Voronoi cell, but you can reconstructed it straightforwardly. * For instance, to find all the edges of Voronoi cell 1, search the output .edge file * for every edge that has input vertex 1 as an endpoint. The corresponding dual * edges in the output .v.edge file form the boundary of Voronoi cell 1." */ uint64_t dim[GMT_DIM_SIZE] = {1, 0, 2, 2}; int i, j, k, n, km1, j2, i2, seg, n_int_edges, n_edges, first_edge = 0, n_extra = 0; int n_to_clip = 0, n_int_vertex = 0, p = 0, corners = 0, n_vertex, change, n_edges_2; unsigned int geometry, side, corner; char header[GMT_LEN64] = {""}; unsigned char *point_type = NULL; struct triangulateio In, Out, vorOut; struct GMT_DATASET *P = NULL; struct GMT_DATASEGMENT *S = NULL; double dy, new_x, xe, ye, xp, yp, x0, y0; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Voronoi partitioning calculated by Jonathan Shewchuk's Triangle [http://www.cs.cmu.edu/~quake/triangle.html]\n"); /* Set everything to 0 and NULL */ gmt_M_memset (&In, 1, struct triangulateio); gmt_M_memset (&Out, 1, struct triangulateio); gmt_M_memset (&vorOut, 1, struct triangulateio); /* Allocate memory for input points */ In.numberofpoints = n = (int)n64; In.pointlist = gmt_M_memory (GMT, NULL, 2 * n, double); /* Copy x,y points to In structure array */ for (i = j = 0; i < n; i++) { In.pointlist[j++] = x_in[i]; In.pointlist[j++] = y_in[i]; } /* Call Jonathan Shewchuk's triangulate algorithm. This is 64-bit safe since * all the structures use 4-byte ints (longs are used internally). The options are * z : Start numbering at zero instead of 1. * I : No iteration numbers * Q : Quiet - suppress all explanations * B : Suppresses output of boundary information * v : Get Voronoi vertices. * D : Conforming Delaunay: all triangles are truly Delaunay. * j : jettison unused vertices from list. */ triangulate ("zIQBvDj", &In, &Out, &vorOut); /* Determine output size for all edges */ n_int_edges = vorOut.numberofedges; /* Count Voronoi vertices and number of infinite rays */ for (i = 0, k = 0; i < n_int_edges; i++, k += 2) { if (vorOut.edgelist[k+1] == -1) n_extra++; /* Infinite rays */ if (vorOut.edgelist[k] > n_int_vertex) n_int_vertex = vorOut.edgelist[k]; if (vorOut.edgelist[k+1] > n_int_vertex) n_int_vertex = vorOut.edgelist[k+1]; } /* Count Voronoi vertices outside w/e/s/n region */ for (i = k = 0; i < n_int_vertex; i++, k += 2) if (vorOut.pointlist[k] < wesn[XLO] || vorOut.pointlist[k] > wesn[XHI] || vorOut.pointlist[k+1] < wesn[YLO] || vorOut.pointlist[k+1] > wesn[YHI]) n_to_clip++; #ifdef DEBUG GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Output from triangulate:\n"); for (i = k = 0; i < n_int_edges; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d normlist = %8g\t%8g\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1], vorOut.normlist[k], vorOut.normlist[k+1]); for (i = k = 0; i < n_int_vertex; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Point %8" PRIu64 " at %g\t%g\n", i, vorOut.pointlist[k], vorOut.pointlist[k+1]); #endif GMT_Report (GMT->parent, GMT_MSG_DEBUG, "n_vertex = %d n_edge = %d n_inf_rays = %d n_outside = %d\n", n_int_vertex, n_int_edges, n_extra, n_to_clip); n_int_vertex++; /* The next edge number */ p = 2 * n_int_vertex; /* Index into vorOut.pointlist for next point to be added along boundary */ n_vertex = n_int_vertex; /* Number of all vertices so far */ if (mode) { /* Want closed Voronoi polygons, so must allocate more space to hold boundary rays and corners */ corners = 4; /* Need to add this many boundingbox corner coordinates */ dim[GMT_SEG] = n; /* When all Voronoi polygons are closed there will be one surrounding each input data point */ dim[GMT_ROW] = 0; /* Variable row length so cannot specify it here - allocate later */ geometry = GMT_IS_POLY; /* Since we will be making closed polygons */ /* Allocate array for point type. This holds which side we exited. Normally 0-3, here we add 1 to use 1-4 * instead since here, 0 means interior point. 1 is south and then we go CCW to 2 (east), 3 (north) and 4 (west) */ point_type = gmt_M_memory (GMT, NULL, n_int_edges + n_extra + corners + n_to_clip, char); } else { /* Want line edges only, all dimensions are already known */ dim[GMT_SEG] = n_int_edges; /* All dimensions are known since we issue just 1 line per segment */ geometry = GMT_IS_LINE; /* Since we only report edges */ } /* Create dataset with a single table with one segment per edge or polygon */ if ((P = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, geometry, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create a data set for gmtsupport_voronoi_shewchuk\n"); gmt_M_free (GMT, point_type); GMT->parent->error = GMT_RUNTIME_ERROR; return NULL; } /* Reallocate the triangle arrays to hold the extra vertices we will need to add */ vorOut.pointlist = realloc (vorOut.pointlist, 2 * (n_vertex + n_extra + corners + 2*n_to_clip) * sizeof (double)); vorOut.edgelist = realloc (vorOut.edgelist, 2 * (n_int_edges + n_extra + corners + n_to_clip) * sizeof (int)); /* First replace infinite rays with finite points where they intersect the border (i.e., we clip the rays to w/e/s/n) */ for (i = k = 0; i < n_int_edges; i++, k++) { km1 = k++; /* Index to this edgelist's first point */ j2 = 2 * vorOut.edgelist[km1]; /* Index into pointlist with x-coordinate of this edge point */ /* Get the coordinates of the point P; this could be an interior [keep] or exterior Voronoi vertex [to be skipped] */ xp = vorOut.pointlist[j2++]; yp = vorOut.pointlist[j2]; if (vorOut.edgelist[k] == -1) { /* Infinite ray; calculate intersection with region boundary */ /* Each edgelist always has a Voronoi vertex as the first point so j2 is never 2 * (-1) */ if (xp < wesn[XLO] || xp > wesn[XHI] || yp < wesn[YLO] || yp > wesn[YHI]) { /* Infinite ray outside boundary - skip edge */ vorOut.edgelist[km1] = -1; /* Mark as a skipped edge */ continue; } /* Determine (xe, ye), the intersection of the ray and the bounding box */ if (vorOut.normlist[km1] < 0.0) { /* Ray will intersect x = xmin, called side = 4 */ xe = wesn[XLO]; side = 4; } else { /* Ray will intersect x = xmax, called side = 2 */ xe = wesn[XHI]; side = 2; } /* Compute y-value at the intersection or ray and border */ dy = fabs ((vorOut.normlist[k] / vorOut.normlist[km1]) * (xe - xp)); ye = yp + dy * copysign (1.0, vorOut.normlist[k]); if (ye < wesn[YLO]) { /* Recompute the x-crossing along y = ymin instead and set ye to ymin */ side = 1; /* South */ new_x = xp + (wesn[YLO] - yp) * (xe - xp) / (ye - yp); GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Voronoi infinite ray truncated from %g %g to %g %g\n", xe, ye, new_x, wesn[YLO]); xe = new_x; ye = wesn[YLO]; } else if (ye > wesn[YHI]) { /* Recompute the x-crossing along y = ymax instead and set ye to ymax */ side = 3; /* North */ new_x = xp + (wesn[YHI] - yp) * (xe - xp) / (ye - yp); GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Voronoi infinite ray truncated from %g %g to %g %g\n", xe, ye, new_x, wesn[YHI]); xe = new_x; ye = wesn[YHI]; } /* Update the truncated ray (-1) in the edge list with a new vertex and add the vertex coordinates to pointlist */ if (mode) point_type[n_vertex] = (unsigned char)side; /* Mark as a border point 1-4 */ vorOut.edgelist[k] = n_vertex++; /* Replace the -1 with the actual point on the boundary */ vorOut.pointlist[p++] = xe; /* Add the ray intersection point to the pointlist */ vorOut.pointlist[p++] = ye; } else { /* A regular edge specified by two points */ i2 = 2 * vorOut.edgelist[k]; /* 2nd index into pointlist with x-coordinate */ /* Get the coordinates of the second Voronoi vertex */ x0 = vorOut.pointlist[i2++]; y0 = vorOut.pointlist[i2]; /* Must check if one of the two points (xp,yp) and (x0,y0) lies outside the region; if so compute intersection with side */ if (xp < wesn[XLO]) { /* Crossing at xmin */ xe = wesn[XLO]; ye = y0 - (y0 - yp) * (x0 - xe) / (x0 - xp); change = k - 1; side = 4; } else if (xp > wesn[XHI]) { /* Crossing at xmax */ xe = wesn[XHI]; ye = y0 - (y0 - yp) * (x0 - xe) / (x0 - xp); change = k - 1; side = 2; } else if (yp < wesn[YLO]) { /* Crossing at ymin */ ye = wesn[YLO]; xe = x0 - (x0 - xp) * (y0 - ye) / (y0 - yp); change = k - 1; side = 1; } else if (yp > wesn[YHI]) { /* Crossing at ymax */ ye = wesn[YHI]; xe = x0 - (x0 - xp) * (y0 - ye) / (y0 - yp); change = k - 1; side = 3; } else if (x0 < wesn[XLO]) { /* Crossing at xmin */ xe = wesn[XLO]; ye = yp - (yp - y0) * (xp - xe) / (xp - x0); change = k; side = 4; } else if (x0 > wesn[XHI]) { /* Crossing at xmax */ xe = wesn[XHI]; ye = yp - (yp - y0) * (xp - xe) / (xp - x0); change = k; side = 2; } else if (y0 < wesn[YLO]) { /* Crossing at ymin */ ye = wesn[YLO]; xe = xp - (xp - x0) * (yp - ye) / (yp - y0); change = k; side = 1; } else if (y0 > wesn[YHI]) { /* Crossing at ymax */ ye = wesn[YHI]; xe = xp - (xp - x0) * (yp - ye) / (yp - y0); change = k; side = 3; } else /* Normal edge - nothing to do for now */ continue; /* Here we replace the edge vertex point with the intersection point and add that point as a new point */ if (mode) point_type[n_vertex] = (unsigned char)side; /* Mark new point as a border point 1-4 */ vorOut.edgelist[change] = n_vertex++; /* Update edgelist with new point on the border, then increase point count */ vorOut.pointlist[p++] = xe; /* Place the new coordinates into the pointlist array */ vorOut.pointlist[p++] = ye; } } /* Now edgelist only contains actual point IDs and pointlist has been updated to hold all added border points */ #ifdef DEBUG GMT_Report (GMT->parent, GMT_MSG_DEBUG, "After infinite ray and exterior vertex crossing replacements:\n"); for (i = k = 0; i < n_int_edges; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d normlist = %8g\t%8g\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1], vorOut.normlist[k], vorOut.normlist[k+1]); for (i = k = 0; i < n_vertex; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Point %8" PRIu64 " at %g\t%g\n", i, vorOut.pointlist[k], vorOut.pointlist[k+1]); #endif if (mode) { /* Need to make closed polygons from edges */ bool first_turn, go_i, go_j; signed char *edge_use = NULL; int expected_sign = 0, prev_sign, next_sign, *start_vertex = NULL, *stop_vertex = NULL; int64_t edge, pstart, pstop, start_point, next_point, prev_point, np; double xstart, xstop, xend, ystart, ystop, yend, cross_product; double *xcoord = NULL, *ycoord = NULL, *dx = NULL, *dy = NULL, *x_pos = NULL, *y_pos = NULL; /* Add the 4 border corners to the point list as well and mark them as corners with special side = 5-8 (i.e., 4 + side 1-4). We need different sides for corners since xmax,ymax is the end for two separate axes. */ vorOut.pointlist[p++] = wesn[XHI]; vorOut.pointlist[p++] = wesn[YLO]; point_type[n_vertex++] = 5; vorOut.pointlist[p++] = wesn[XHI]; vorOut.pointlist[p++] = wesn[YHI]; point_type[n_vertex++] = 6; vorOut.pointlist[p++] = wesn[XLO]; vorOut.pointlist[p++] = wesn[YHI]; point_type[n_vertex++] = 7; vorOut.pointlist[p++] = wesn[XLO]; vorOut.pointlist[p++] = wesn[YLO]; point_type[n_vertex++] = 8; #ifdef DEBUG GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Before border edges are added:\n"); for (i = k = 0; i < n_int_edges; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d normlist = %8g\t%8g\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1], vorOut.normlist[k], vorOut.normlist[k+1]); for (i = k = 0; i < n_vertex; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Point %8" PRIu64 " [%d] at %g\t%g\n", i, point_type[i], vorOut.pointlist[k], vorOut.pointlist[k+1]); #endif /* Finally add the new edges between new points along the border and with the corners */ /* Do the y = ymin and y = ymax sides first (i.e., side 1 and 3) */ edge = 2 * n_int_edges; /* Start index of next added edge in edgelist */ for (side = 1; side <= 3; side += 2) { /* Since almost the same logic for the two sides */ corner = (side == 1) ? 5 : 6; /* The corner at the end of this x-axis */ xstart = wesn[XLO]; pstart = (side == 1) ? n_vertex - 1 : n_vertex - 2; /* Initialize the start of the edge at beginning of x-axis */ xstop = xend = wesn[XHI]; pstop = (side == 1) ? n_vertex - 4 : n_vertex - 3; /* End of x-axis is end of edge unless there are intersections */ while (xstart < wesn[XHI]) { /* As long as we have incomplete edges */ for (i = n_int_vertex; i < n_vertex; i++) { /* Must check all of them for one that could be a candidate along this x-axis */ if (!(point_type[i] == side || point_type[i] == corner)) continue; /* Not along the x-axis including the end corner */ j2 = 2 * i; /* Index into pointlist for the x-coordinate */ if (vorOut.pointlist[j2] <= xstop && vorOut.pointlist[j2] > xstart) { /* An intersection closer to the start of our edge, update xstop */ xstop = vorOut.pointlist[j2]; pstop = i; } } /* Here, xstop is the point closest to xstart and pstop is its index. */ /* Add this edge to edgelist */ vorOut.edgelist[edge++] = (int)pstart; vorOut.edgelist[edge++] = (int)pstop; pstart = pstop; xstart = xstop; xstop = xend; /* Let the end of this edge become start of the next edge */ } } /* Do the xmin and xmax sides next (sides 2 and 4) */ for (side = 2; side <= 4; side += 2) { corner = (side == 2) ? 6 : 7; /* The corner at the end of this y-axis */ ystart = wesn[YLO]; pstart = (side == 2) ? n_vertex - 4 : n_vertex - 1; /* Initialize the start of the edge at beginning of y-axis */ ystop = yend = wesn[YHI]; pstop = (side == 2) ? n_vertex - 3 : n_vertex - 2; /* End of y-axis is end of edge unless there are intersections */ while (ystart < wesn[YHI]) { /* As long as we have incomplete edges */ for (i = n_int_vertex; i < n_vertex; i++) { /* Must check all of them for one that could be a candidate along this y-axis */ if (!(point_type[i] == side || point_type[i] == corner)) continue; /* Not along the y-axis including the end corner */ j2 = 2 * i + 1; /* Index into pointlist for the y-coordinate */ if (vorOut.pointlist[j2] <= ystop && vorOut.pointlist[j2] > ystart) { /* An intersection closer to the start of our edge, update ystop */ ystop = vorOut.pointlist[j2]; pstop = i; } } /* Here, ystop is the point closest to xstart and pstop is its index. */ /* Add this edge to edgelist */ vorOut.edgelist[edge++] = (int)pstart; vorOut.edgelist[edge++] = (int)pstop; pstart = pstop; ystart = ystop; ystop = yend; /* Let the end of this edge become start of the next edge */ } } n_edges = (int)edge; /* Total number of all edges times 2 */ #ifdef DEBUG GMT_Report (GMT->parent, GMT_MSG_DEBUG, "\nAfter border edges are added:\n"); for (i = k = 0; k < n_edges; i++, k += 2) { if (i < n_int_edges) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d normlist = %8g\t%8g\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1], vorOut.normlist[k], vorOut.normlist[k+1]); else GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1]); } for (i = k = 0; i < n_vertex; i++, k += 2) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Point %8" PRIu64 " [%d] at %g\t%g\n", i, point_type[i], vorOut.pointlist[k], vorOut.pointlist[k+1]); #endif gmt_M_free (GMT, point_type); vorOut.normlist = realloc (vorOut.normlist, n_edges * sizeof (double)); /* Remove the unneeded external edges flagged by two -1 signs in the edgelist */ for (j2 = i2 = 0; j2 < n_edges; j2 += 2) { /* For all the edges */ if (i2 < j2) { /* Shuffle the entries to lower indices */ vorOut.edgelist[i2] = vorOut.edgelist[j2]; vorOut.edgelist[i2+1] = vorOut.edgelist[j2+1]; vorOut.normlist[i2] = vorOut.normlist[j2]; vorOut.normlist[i2+1] = vorOut.normlist[j2+1]; } if (!(vorOut.edgelist[j2] == -1 && vorOut.edgelist[j2+1] == -1)) i2 += 2; /* Only increment output index when not a dummy edge */ } /* Update the count of edges. All the unneeded edges will lie within the original n_int_edges section */ n_edges -= (j2 - i2); n_int_edges -= ((j2 - i2)/2); #ifdef DEBUG GMT_Report (GMT->parent, GMT_MSG_DEBUG, "\nAfter removing unused edges:\n"); for (i = k = 0; k < n_edges; i++, k += 2) { if (i < n_int_edges) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d normlist = %8g\t%8g\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1], vorOut.normlist[k], vorOut.normlist[k+1]); else GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Edge %8" PRIu64 " Point %8d to %8d\n", i, vorOut.edgelist[k], vorOut.edgelist[k+1]); } #endif /* Precalculate dx, dy for each edge and store these in normlist (which we need to reallocate first) */ n_edges_2 = n_edges / 2; edge_use = gmt_M_memory (GMT, NULL, n_edges_2, signed char); /* 0 = unused, +/-1 = used once in normal or reverse direction, 2 = used twice and done with */ dx = vorOut.normlist; dy = &vorOut.normlist[1]; /* So we can use dx[index] and dy[index] */ start_vertex = vorOut.edgelist; stop_vertex = &vorOut.edgelist[1]; x_pos = vorOut.pointlist; y_pos = &vorOut.pointlist[1]; for (i = k = 0; k < n_edges; i++, k += 2) { /* For all the edges */ i2 = 2 * start_vertex[k]; /* Location of x and y coordinates for start of edge */ j2 = 2 * stop_vertex[k]; /* Location of x and y coordinates for end of edge */ dx[k] = x_pos[j2] - x_pos[i2]; /* dx = end_x - start_x */ dy[k] = y_pos[j2] - y_pos[i2]; /* dy = end_y - start_y */ if (i >= n_int_edges) edge_use[i] = 1; /* These edges will only be used once (sign is not important here) */ } /* Need temp array to hold the coordinates of a single polygon [I hope this is large enough - don't know yet] */ xcoord = gmt_M_memory (GMT, NULL, n_int_edges/2, double); ycoord = gmt_M_memory (GMT, NULL, n_int_edges/2, double); /* Now stitch together the closed polygons */ seg = 0; first_turn = true; /* When we find the first 3 points in a polygon we can compute the cross-product and determine sign for all polygons */ while (seg < n) { /* Find the first interior edge that has not already been used twice */ for (i = first_edge, go_i = true; go_i && i < n_int_edges; i++) { if (edge_use[i] == 2) { /* Done with this edge */ first_edge++; /* Might as well skip this one from now on */ continue; } /* Here, i2 = 2*i is the index into the edgelist and normlist arrays for the previous edge */ i2 = 2 * i; if (edge_use[i] == +1) { /* Must reverse direction the 2nd time we start with this edge */ start_point = prev_point = stop_vertex[i2]; next_point = start_vertex[i2]; prev_sign = -1; /* Since we are reversing this edge */ edge_use[i] = 2; /* Flag we are finished with this edge */ } else if (edge_use[i] == -1) { /* Go in the forward direction the 2nd time we start with this edge */ start_point = prev_point = start_vertex[i2]; next_point = stop_vertex[i2]; prev_sign = +1; /* Since we are NOT reversing this edge */ edge_use[i] = 2; /* Flag we are finished with this edge */ } else { /* Go in the forward direction the first time */ start_point = prev_point = start_vertex[i2]; next_point = stop_vertex[i2]; prev_sign = +1; /* Since we are NOT reversing this edge */ edge_use[i] = 1; /* Flag we have used this edge one more time */ } /* Add this edge as the first line in the new polygon */ xcoord[0] = x_pos[2*start_point]; ycoord[0] = y_pos[2*start_point]; xcoord[1] = x_pos[2*next_point]; ycoord[1] = y_pos[2*next_point]; np = 2; /* Total points in this polygon so far */ while (next_point != start_point) { /* When this is false we have automatically closed the polygon */ /* Find first edge of all edges not equal the current edge that has next_point either as 1st or 2nd vertex and the other vertex is not our prev_point */ for (j = first_edge, go_j = true; go_j && j < n_edges_2; j++) { if (j == i || edge_use[j] == 2) continue; /* Either the same edge or we are already done with this edge */ /* Here, j2 = 2*j is the index into the edgelist and normlist arrays for the next edge */ j2 = 2 * j; if (start_vertex[j2] == next_point && stop_vertex[j2] != prev_point) /* Must go in forward order */ next_sign = 1; else if (stop_vertex[j2] == next_point && start_vertex[j2] != prev_point) /* Go in reverse order */ next_sign = -1; else /* Not a desired edge since we could not match the vertex number */ continue; /* Here we found an edge candidate, but we don't know if it leads to a convex angle */ cross_product = dx[i2] * dy[j2] - dx[j2] * dy[i2]; /* This should be negative for right turns */ if (first_turn) { /* First pair of edges defines the desired sign of all the subsequent cross-products */ expected_sign = prev_sign * next_sign * gmt_signum (cross_product); /* sign may flip if any edge was reversed */ first_turn = false; } else if ((prev_sign * next_sign * gmt_signum (cross_product)) != expected_sign) /* Not making a convex polygon */ continue; /* Here we are going in a convex direction, so we add this new edge to the polygon */ prev_point = next_point; next_point = (next_sign == 1) ? stop_vertex[j2] : start_vertex[j2]; prev_sign = next_sign; i2 = j2; /* Add next_point to our polygon */ xcoord[np] = x_pos[2*next_point]; ycoord[np] = y_pos[2*next_point]; if (edge_use[j]) /* Now used twice */ edge_use[j] = 2; else /* Flag according to which direction we used it */ edge_use[j] = (char)next_sign; np++; /* Increase polygon length counter */ go_j = false; /* Exit this loop and search for next edge */ } } /* Finalize the polygon */ snprintf (header, GMT_LEN64, "Voronoi polygon # %d -L%d", seg, seg); S = P->table[0]->segment[seg] = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, np, 2U, header, P->table[0]->segment[seg]); gmt_M_memcpy (S->data[GMT_X], xcoord, np, double); gmt_M_memcpy (S->data[GMT_Y], ycoord, np, double); P->table[0]->n_records += np; /* Update counts */ P->n_records += np; seg++; go_i = false; /* Exit loop and move to the next polygon */ } } /* Free temporary arrays */ gmt_M_free (GMT, xcoord); gmt_M_free (GMT, ycoord); gmt_M_free (GMT, edge_use); } else { /* Only want line edges */ for (seg = k = 0; seg < n_int_edges; seg++) { j2 = 2 * vorOut.edgelist[k++]; S = P->table[0]->segment[seg]; /* Current output edge segment */ S->data[GMT_X][0] = vorOut.pointlist[j2++]; S->data[GMT_Y][0] = vorOut.pointlist[j2]; j2 = 2 * vorOut.edgelist[k++]; S->data[GMT_X][1] = vorOut.pointlist[j2++]; S->data[GMT_Y][1] = vorOut.pointlist[j2]; snprintf (header, GMT_LEN64, "Voronoi edge # %d -L%d", seg, seg); S->header = strdup (header); } } gmt_set_dataset_minmax (GMT, P); /* Determine min/max for each column */ /* Free the triangulate arrays that were all allocated internally */ gmt_M_str_free (Out.pointlist); gmt_M_str_free (Out.trianglelist); gmt_M_str_free (vorOut.pointattributelist); gmt_M_str_free (vorOut.pointlist); gmt_M_str_free (vorOut.edgelist); gmt_M_str_free (vorOut.normlist); gmt_M_free (GMT, In.pointlist); return (P); } #else /*! Dummy functions since not installed */ GMT_LOCAL uint64_t gmtsupport_delaunay_shewchuk (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n, int **link) { gmt_M_unused (x_in); gmt_M_unused (y_in); gmt_M_unused (n); gmt_M_unused (link); GMT_Report (GMT->parent, GMT_MSG_ERROR, "unavailable: Shewchuk's triangle option was not selected during GMT installation\n"); return (0); } /*! . */ GMT_LOCAL struct GMT_DATASET * gmtsupport_voronoi_shewchuk (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n, double *wesn, unsigned int mode) { gmt_M_unused (x_in); gmt_M_unused (y_in); gmt_M_unused (n); gmt_M_unused (wesn); gmt_M_unused (mode); GMT_Report (GMT->parent, GMT_MSG_ERROR, "unavailable: Shewchuk's triangle option was not selected during GMT installation\n"); return (NULL); } #endif /* * gmt_delaunay performs a Delaunay triangulation on the input data * and returns a list of indices of the points for each triangle * found. Algorithm translated from * Watson, D. F., ACORD: Automatic contouring of raw data, * Computers & Geosciences, 8, 97-101, 1982. */ struct GMT_PAIR { double x, y; uint64_t rec; }; /*! . */ GMT_LOCAL int gmtsupport_sort_pair (const void *p_1, const void *p_2) { const struct GMT_PAIR *point_1 = (const struct GMT_PAIR *)p_1, *point_2 = (const struct GMT_PAIR *)p_2; if (point_1->x < point_2->x) return -1; if (point_1->x > point_2->x) return +1; if (point_1->y < point_2->y) return -1; if (point_1->y > point_2->y) return +1; return 0; } /* Leave link as int**, not int** */ /*! . */ GMT_LOCAL uint64_t gmtsupport_delaunay_watson (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n, int **link) { /* Input point x coordinates */ /* Input point y coordinates */ /* Number of input points */ /* pointer to List of point ids per triangle. Vertices for triangle no i is in link[i*3], link[i*3+1], link[i*3+2] */ int *index = NULL; /* Must be int not uint64_t */ int ix[3], iy[3]; bool done; uint64_t i, j, nuc, ij, jt, km, id, isp, l1, l2, k, k1, jz, i2, kmt, kt, size; int64_t *istack = NULL, *x_tmp = NULL, *y_tmp = NULL; double det[2][3], *x_circum = NULL, *y_circum = NULL, *r2_circum = NULL, *x = NULL, *y = NULL; double xmin, xmax, ymin, ymax, datax, dx, dy, dsq, dd; GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Delaunay triangulation calculated by Dave Watson's ACORD [Computers & Geosciences, 8, 97-101, 1982]\n"); { /* Note 2019/01/07: We were notified via https://github.com/GenericMappingTools/gmt/issues/279 * that the Watson algorithm may give junk if there are duplicate entries in the input, and if so we issue * a stern warning to users so they can clean up the file first before calling gmtsupport_delaunay_watson */ struct GMT_PAIR *P = gmt_M_memory (GMT, NULL, n, struct GMT_PAIR); uint64_t n_duplicates = 0; for (i = 0; i < n; i++) { P[i].x = x_in[i]; P[i].y = y_in[i]; P[i].rec = i + 1; } qsort (P, n, sizeof (struct GMT_PAIR), gmtsupport_sort_pair); for (i = 1; i < n; i++) { if (doubleAlmostEqualZero (P[i].x, P[i-1].x) && doubleAlmostEqualZero (P[i].y, P[i-1].y)) { GMT_Report (GMT->parent, GMT_MSG_WARNING, "Records %" PRIu64 " and %" PRIu64 " are duplicates!\n", P[i-1].rec, P[i].rec); n_duplicates++; } } if (n_duplicates) { /* Clearly not monotonically increasing or decreasing in x */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Bug Report Advice for Watson ACORD External Code:\n"); GMT_Report (GMT->parent, GMT_MSG_WARNING, "The Watson algorithm can fail if there are duplicate (x,y) records.\n"); GMT_Report (GMT->parent, GMT_MSG_WARNING, "We found %" PRIu64 " duplicate records in your data set.\n", n_duplicates); GMT_Report (GMT->parent, GMT_MSG_WARNING, "Please remove duplicates and redo your analysis if the results are corrupted.\n"); } gmt_M_free (GMT, P); } size = 10 * n + 1; n += 3; index = gmt_M_memory (GMT, NULL, 3 * size, int); istack = gmt_M_memory (GMT, NULL, size, int64_t); x_tmp = gmt_M_memory (GMT, NULL, size, int64_t); y_tmp = gmt_M_memory (GMT, NULL, size, int64_t); x_circum = gmt_M_memory (GMT, NULL, size, double); y_circum = gmt_M_memory (GMT, NULL, size, double); r2_circum = gmt_M_memory (GMT, NULL, size, double); x = gmt_M_memory (GMT, NULL, n, double); y = gmt_M_memory (GMT, NULL, n, double); x[0] = x[1] = -1.0; x[2] = 5.0; y[0] = y[2] = -1.0; y[1] = 5.0; x_circum[0] = y_circum[0] = 2.0; r2_circum[0] = 18.0; ix[0] = ix[1] = 0; ix[2] = 1; iy[0] = 1; iy[1] = iy[2] = 2; for (i = 0; i < 3; i++) index[i] = (int)i; for (i = 0; i < size; i++) istack[i] = i; xmin = ymin = 1.0e100; xmax = ymax = -1.0e100; for (i = 3, j = 0; i < n; i++, j++) { /* Copy data and get extremas */ x[i] = x_in[j]; y[i] = y_in[j]; if (x[i] > xmax) xmax = x[i]; if (x[i] < xmin) xmin = x[i]; if (y[i] > ymax) ymax = y[i]; if (y[i] < ymin) ymin = y[i]; } /* Normalize data */ datax = 1.0 / MAX (xmax - xmin, ymax - ymin); for (i = 3; i < n; i++) { x[i] = (x[i] - xmin) * datax; y[i] = (y[i] - ymin) * datax; } isp = id = 1; for (nuc = 3; nuc < n; nuc++) { km = 0; for (jt = 0; jt < isp; jt++) { /* loop through established 3-tuples */ ij = 3 * jt; /* Test if new data is within the jt circumcircle */ dx = x[nuc] - x_circum[jt]; if ((dsq = r2_circum[jt] - dx * dx) < 0.0) continue; dy = y[nuc] - y_circum[jt]; if ((dsq -= dy * dy) < 0.0) continue; /* Delete this 3-tuple but save its edges */ id--; istack[id] = jt; /* Add edges to x/y_tmp but delete if already listed */ for (i = 0; i < 3; i++) { l1 = ix[i]; l2 = iy[i]; if (km > 0) { kmt = km; for (j = 0, done = false; !done && j < kmt; j++) { if (index[ij+l1] != x_tmp[j]) continue; if (index[ij+l2] != y_tmp[j]) continue; km--; if (j >= km) { done = true; continue; } for (k = j; k < km; k++) { k1 = k + 1; x_tmp[k] = x_tmp[k1]; y_tmp[k] = y_tmp[k1]; done = true; } } } else done = false; if (!done) { x_tmp[km] = index[ij+l1]; y_tmp[km] = index[ij+l2]; km++; } } } /* Form new 3-tuples */ for (i = 0; i < km; i++) { kt = istack[id]; ij = 3 * kt; id++; /* Calculate the circumcircle center and radius * squared of points ktetr[i,*] and place in tetr[kt,*] */ for (jz = 0; jz < 2; jz++) { i2 = (jz == 0) ? x_tmp[i] : y_tmp[i]; det[jz][0] = x[i2] - x[nuc]; det[jz][1] = y[i2] - y[nuc]; det[jz][2] = 0.5 * (det[jz][0] * (x[i2] + x[nuc]) + det[jz][1] * (y[i2] + y[nuc])); } dd = 1.0 / (det[0][0] * det[1][1] - det[0][1] * det[1][0]); x_circum[kt] = (det[0][2] * det[1][1] - det[1][2] * det[0][1]) * dd; y_circum[kt] = (det[0][0] * det[1][2] - det[1][0] * det[0][2]) * dd; dx = x[nuc] - x_circum[kt]; dy = y[nuc] - y_circum[kt]; r2_circum[kt] = dx * dx + dy * dy; index[ij++] = (int)x_tmp[i]; index[ij++] = (int)y_tmp[i]; index[ij] = (int)nuc; } isp += 2; } for (jt = i = 0; jt < isp; jt++) { ij = 3 * jt; if (index[ij] < 3 || r2_circum[jt] > 1.0) continue; index[i++] = index[ij++] - 3; index[i++] = index[ij++] - 3; index[i++] = index[ij++] - 3; } index = gmt_M_memory (GMT, index, i, int); *link = index; gmt_M_free (GMT, istack); gmt_M_free (GMT, x_tmp); gmt_M_free (GMT, y_tmp); gmt_M_free (GMT, x_circum); gmt_M_free (GMT, y_circum); gmt_M_free (GMT, r2_circum); gmt_M_free (GMT, x); gmt_M_free (GMT, y); return (i/3); } /*! . */ GMT_LOCAL struct GMT_DATASET * gmtsupport_voronoi_watson (struct GMT_CTRL *GMT, double *x_in, double *y_in, uint64_t n, double *wesn, unsigned int mode) { gmt_M_unused(x_in); gmt_M_unused(y_in); gmt_M_unused(n); gmt_M_unused(wesn); gmt_M_unused(mode); GMT_Report (GMT->parent, GMT_MSG_ERROR, "No Voronoi unless you select Shewchuk's triangle option during GMT installation\n"); return (0); } /*! . */ GMT_LOCAL int gmtsupport_ensure_new_mapinset_syntax (struct GMT_CTRL *GMT, char option, char *in_text, char *text, char *panel_txt) { /* Recasts any old syntax using new syntax and gives a warning. Assumes text and panel_text are blank and have adequate space */ if (strstr (in_text, "+c") || strstr (in_text, "+g") || strstr (in_text, "+p")) { /* Tell-tale sign of old syntax */ char p[GMT_BUFSIZ] = {""}, txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}; unsigned int pos = 0, start = 0, i; int n; bool center = false; for (i = 0; start == 0 && in_text[i]; i++) { /* Find the first valid modifier */ if (in_text[i] == '+') { if (strchr ("cgp", in_text[i+1])) start = i; /* Found start of the modifiers */ } } while ((gmt_strtok (&in_text[start], "+", &pos, p))) { switch (p[0]) { case 'c': /* Got inset center */ center = true; if ((n = sscanf (&p[1], "%[^/]/%s", txt_a, txt_b)) != 2) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Must specify +c/ for center [Also note this is obsolete syntax]\n", option); return (1); } break; case 'g': /* Fill specification [Obsolete, now done via panel attributes ] */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Option -%c: Insert fill attributes now given with panel setting (-F)\n", option); strcat (panel_txt, "+g"); strcat (panel_txt, &p[1]); break; case 'p': /* Pen specification [Obsolete, now done via panel attributes ] */ GMT_Report (GMT->parent, GMT_MSG_WARNING, "Option -%c: Insert pen attributes now given with panel setting (-F)\n", option); strcat (panel_txt, "+p"); strcat (panel_txt, &p[1]); break; default: break; } } in_text[start] = '\0'; /* Chop off modifiers for now */ if (center) { /* Must extract dimensions of map inset */ char unit[2] = {0, 0}; sprintf (text, "g%s/%s/", txt_a, txt_b); /* -Dg/ is the new reference point */ n = sscanf (in_text, "%[^/]/%s", txt_a, txt_b); /* Read dimensions */ if (strchr (GMT_LEN_UNITS2, txt_a[0])) { /* Dimensions in allowable distance units, with unit in front of distance */ unit[0] = txt_a[0]; /* Extract the unit */ strcat (text, "+w"); /* Append width modifier */ strcat (text, &txt_a[1]); /* Append width to new option */ strcat (text, unit); /* Append unit to new option */ if (n == 2) { /* Got separate height */ strcat (text, "/"); /* Separate width from height */ strcat (text, txt_b); /* Append height or duplicate */ strcat (text, unit); /* Append unit */ } } else strcat (text, in_text); /* Append h/w as is */ strcat (text, "+jCM"); /* Append justification */ } else /* Just keep as is, minus the modifiers */ strcpy (text, in_text); in_text[start] = '+'; /* Restore modifiers */ GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Converted %s to %s and %s\n", in_text, text, panel_txt); } else /* New syntax or no +g,+p given so old syntax is same as new */ strcpy (text, in_text); return (0); } /*! . */ GMT_LOCAL int gmtsupport_getscale_old (struct GMT_CTRL *GMT, char option, char *text, struct GMT_MAP_SCALE *ms) { /* This function parses the -L map scale syntax: * -L[f][x]/[/]//[e|f|M|n|k|u][+u] * The function is also backwards compatible with the previous map scale syntax: * -L [f][x]/[/]//[e|f|M|n|k|u][+l