1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Author: Paul Wessel
19 * Date: 1-JAN-2010
20 * Version: 6 API
21 *
22 * Brief synopsis: grdimage will read a grid or image and plot the area
23 * with optional shading using the selected projection. Optionally, it
24 * may return the raw raster image.
25 *
26 */
27
28 #include "gmt_dev.h"
29
30 #define THIS_MODULE_CLASSIC_NAME "grdimage"
31 #define THIS_MODULE_MODERN_NAME "grdimage"
32 #define THIS_MODULE_LIB "core"
33 #define THIS_MODULE_PURPOSE "Project and plot grids or images"
34 #define THIS_MODULE_KEYS "<G{+,CC(,IG(,>X},>IA,<ID"
35 #define THIS_MODULE_NEEDS "Jg"
36 #define THIS_MODULE_OPTIONS "->BJKOPRUVXYfnptxy" GMT_OPT("Sc") GMT_ADD_x_OPT
37
38 #ifdef HAVE_GDAL
39 /* These are images that GDAL knows how to read for us. Otherwise we can only deal with ppm */
40 #define N_IMG_EXTENSIONS 6
41 static char *gdal_ext[N_IMG_EXTENSIONS] = {"tiff", "tif", "gif", "png", "jpg", "bmp"};
42 #endif
43
44 #define GRDIMAGE_NAN_INDEX (GMT_NAN - 3)
45
46 /* Control structure for grdimage */
47
48 struct GRDIMAGE_CTRL {
49 struct GRDIMAGE_In { /* Input grid or image */
50 bool active;
51 char *file;
52 } In;
53 struct GRDIMAGE_Out { /* Output image under -A */
54 bool active;
55 char *file;
56 } Out;
57 struct GRDIMAGE_C { /* -C<cpt> or -C<color1>,<color2>[,<color3>,...][+i<dz>] */
58 bool active;
59 double dz; /* Rounding for min/max determined from data */
60 char *file;
61 } C;
62 struct GRDIMAGE_D { /* -D to read image instead of grid */
63 bool active;
64 bool mode; /* Use info of -R option to reference image */
65 } D;
66 struct GRDIMAGE_A { /* -A to write a raster file or return image to API */
67 bool active;
68 bool return_image;
69 unsigned int way;
70 char *file;
71 } A;
72 struct GRDIMAGE_E { /* -Ei|<dpi> */
73 bool active;
74 bool device_dpi;
75 unsigned int dpi;
76 } E;
77 struct GRDIMAGE_G { /* -G<rgb>[+b|f] */
78 bool active;
79 double rgb[2][4];
80 } G;
81 struct GRDIMAGE_I { /* -I[<intens_or_zfile>|<value>|<modifiers>] */
82 bool active;
83 bool constant;
84 bool derive;
85 double value;
86 char *azimuth; /* Default azimuth(s) for shading */
87 char *file;
88 char *method; /* Default scaling method */
89 char *ambient; /* Default ambient offset */
90 } I;
91 struct GRDIMAGE_M { /* -M */
92 bool active;
93 } M;
94 struct GRDIMAGE_N { /* -N */
95 bool active;
96 } N;
97 struct GRDIMAGE_Q { /* -Q[r/g/b] */
98 bool active;
99 bool transp_color; /* true if a color was given */
100 double rgb[4]; /* Pixel value for transparency in images */
101 } Q;
102 struct GRDIMAGE_W { /* -W */
103 bool active;
104 } W;
105 };
106
107 struct GRDIMAGE_CONF {
108 /* Configuration structure for things to pass around to sub-functions */
109 unsigned int colormask_offset; /* Either 0 or 3 depending on -Q */
110 bool int_mode; /* true if we are to apply illumination */
111 unsigned int *actual_row; /* Array with actual rows as a function of pseudo row */
112 unsigned int *actual_col; /* Array of actual columns as a function of pseudo col */
113 int64_t n_columns, n_rows; /* Signed dimensions for use in OpenMP loops */
114 uint64_t nm; /* Number of pixels in the image */
115 struct GMT_PALETTE *P; /* Pointer to the active palette [NULL if image] */
116 struct GMT_GRID *Grid; /* Pointer to the active grid [NULL if image] */
117 struct GMT_GRID *Intens; /* Pointer to the active intensity grid [NULL if no intensity] */
118 struct GMT_IMAGE *Image; /* Pointer to the active image [NUYLL if grid] */
119 };
120
New_Ctrl(struct GMT_CTRL * GMT)121 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
122 struct GRDIMAGE_CTRL *C;
123
124 C = gmt_M_memory (GMT, NULL, 1, struct GRDIMAGE_CTRL);
125
126 /* Initialize values whose defaults are not 0/false/NULL */
127
128 C->G.rgb[GMT_BGD][0] = C->G.rgb[GMT_BGD][1] = C->G.rgb[GMT_BGD][2] = 1.0; /* White background pixel, black foreground defaults */
129 C->I.azimuth = strdup ("-45.0"); /* Default azimuth for shading when -I+d is used */
130 C->I.method = strdup ("t1"); /* Default normalization for shading when -I+d is used */
131 C->I.ambient = strdup ("0"); /* Default ambient light for shading when -I+d is used */
132
133 return (C);
134 }
135
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * C)136 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *C) { /* Deallocate control structure */
137 if (!C) return;
138 gmt_M_str_free (C->In.file);
139 gmt_M_str_free (C->A.file);
140 gmt_M_str_free (C->C.file);
141 gmt_M_str_free (C->I.file);
142 gmt_M_str_free (C->I.azimuth);
143 gmt_M_str_free (C->I.ambient);
144 gmt_M_str_free (C->I.method);
145 gmt_M_str_free (C->Out.file);
146 gmt_M_free (GMT, C);
147 }
148
usage(struct GMTAPI_CTRL * API,int level)149 static int usage (struct GMTAPI_CTRL *API, int level) {
150 #ifdef HAVE_GDAL
151 const char *A = " [-A<out_img>[=<driver>]]";
152 #else
153 const char *A = " [-A]";
154 #endif
155 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
156 const char *extra[2] = {A, " [-A]"};
157 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
158 GMT_Usage (API, 0, "usage: %s %s %s%s [%s] [-C<cpt>] [-D[r]] [-Ei|<dpi>] "
159 "[-G<rgb>[+b|f]] [-I[<intensgrid>|<value>|<modifiers>]] %s[-M] [-N] %s%s[-Q[<color>]] "
160 "[%s] [%s] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s] [%s] [%s]\n",
161 name, GMT_INGRID, GMT_J_OPT, extra[API->external], GMT_B_OPT, API->K_OPT, API->O_OPT, API->P_OPT, GMT_Rgeo_OPT, GMT_U_OPT,
162 GMT_V_OPT, GMT_X_OPT, GMT_Y_OPT, API->c_OPT, GMT_f_OPT, GMT_n_OPT, GMT_p_OPT, GMT_t_OPT, GMT_x_OPT, GMT_PAR_OPT);
163
164 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
165
166 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
167 gmt_ingrid_syntax (API, 0, "Name of a grid or image to plot");
168 GMT_Usage (API, -2, "Note: Grid z-values are in user units and will be "
169 "converted to rgb colors via the CPT. Alternatively, give a raster image. "
170 "If the image is plain (e.g., JPG, PNG, GIF) you must also give a corresponding -R.");
171 if (API->external) /* External interface */
172 GMT_Usage (API, -2, "If -D is used then <grid> is instead expected to be an image.");
173 GMT_Option (API, "J-");
174 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
175 if (API->external) /* External interface */
176 GMT_Usage (API, 1, "\n-A Return a GMT raster image instead of a PostScript plot.");
177 else {
178 #ifdef HAVE_GDAL
179 GMT_Usage (API, 1, "\n-A<out_img>[=<driver>]");
180 GMT_Usage (API, -2, "Save image in a raster format (.bmp, .gif, .jpg, .png, .tif) instead of PostScript. "
181 "If filename does not have any of these extensions then append =<driver> to select "
182 "the desired image format. The 'driver' is the driver code name used by GDAL. "
183 "See GDAL documentation for available drivers. Note: any vector elements are lost.");
184 #else
185 GMT_Usage (API, 1, "\n-A Save image in a PPM format (give .ppm extension) instead of PostScript.");
186 #endif
187 }
188 GMT_Option (API, "B-");
189 GMT_Usage (API, 1, "\n-C<cpt>");
190 GMT_Usage (API, -2, "Color palette file to convert grid values to colors. Optionally, name a master cpt "
191 "to automatically assign continuous colors over the data range [%s]; if so, "
192 "optionally append +i<dz> to quantize the range [the exact grid range]. "
193 "Another option is to specify -C<color1>,<color2>[,<color3>,...] to build a "
194 "linear continuous cpt from those colors automatically.", API->GMT->current.setting.cpt);
195 GMT_Usage (API, 1, "\n-D[r]");
196 if (API->external) /* External interface */
197 GMT_Usage (API, -2, "Input is an image instead of a grid. Append r to equate image region to -R region.");
198 #ifdef HAVE_GDAL
199 else {
200 GMT_Usage (API, -2, "GMT automatically detects standard image formats. For non-standard formats you must "
201 "use -D to force it to be seen as an image. Append r to equate image region with -R region.");
202 }
203 #endif
204 GMT_Usage (API, 1, "\n-Ei|<dpi>");
205 GMT_Usage (API, -2, "Set dpi for the projected grid which must be constructed [100] "
206 "if -J implies a nonlinear graticule [Default gives same size as input grid]. "
207 "Alternatively, append i to do the interpolation in PostScript at device resolution (invalid with -Q).");
208 gmt_rgb_syntax (API->GMT, 'G', "Set transparency color for images that otherwise would result in 1-bit images. ");
209 GMT_Usage (API, 3, "+b Set background color.");
210 GMT_Usage (API, 3, "+f Set foreground color [Default].");
211 GMT_Usage (API, 1, "\n-I[<intensgrid>|<value>|<modifiers>]");
212 GMT_Usage (API, -2, "Apply directional illumination. Append name of an intensity grid, or "
213 "for a constant intensity (i.e., change the ambient light), just give a scalar. "
214 "To derive intensities from <grid> instead, append desired modifiers:");
215 GMT_Usage (API, 3, "+a Specify <azim> of illumination [-45].");
216 GMT_Usage (API, 3, "+n Set the <method> and <scale> to use [t1].");
217 GMT_Usage (API, 3, "+m Set <ambient> light to add [0].");
218 GMT_Usage (API, -2, "Alternatively, use -I+d to accept the default values (see grdgradient for more details). "
219 "To derive intensities from another grid than <grid>, give the alternative data grid with suitable modifiers.");
220 GMT_Option (API, "K");
221 GMT_Usage (API, 1, "\n-M Force a monochrome (gray-scale) image.");
222 GMT_Usage (API, 1, "\n-N Do Not clip image at the map boundary.");
223 GMT_Option (API, "O,P");
224 GMT_Usage (API, 1, "\n-Q[<color>]");
225 GMT_Usage (API, -2, "Use color-masking to make grid nodes with z = NaN or black image pixels transparent. "
226 "Append an alternate <color> to change the transparent pixel for images [black].");
227 GMT_Option (API, "R");
228 GMT_Option (API, "U,V,X,c,f,n,p,t,x,.");
229
230 return (GMT_MODULE_USAGE);
231 }
232
233 /* A few non-exported library functions we need here only */
234
235 EXTERN_MSC int gmtinit_parse_n_option (struct GMT_CTRL *GMT, char *item);
236 EXTERN_MSC int gmtlib_get_grdtype (struct GMT_CTRL *GMT, unsigned int direction, struct GMT_GRID_HEADER *h);
237 EXTERN_MSC int gmtlib_read_grd_info (struct GMT_CTRL *GMT, char *file, struct GMT_GRID_HEADER *header);
238
parse(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GMT_OPTION * options)239 static int parse (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GMT_OPTION *options) {
240 /* This parses the options provided to grdimage and sets parameters in Ctrl.
241 * Note Ctrl has already been initialized and non-zero default values set.
242 * Any GMT common options will override values set previously by other commands.
243 * It also replaces any file names specified as input or output with the data ID
244 * returned when registering these sources/destinations with the API.
245 */
246
247 unsigned int n_errors = 0, n_files = 0, ind, off, k;
248 char *c = NULL, *file[3] = {NULL, NULL, NULL};
249 struct GMT_OPTION *opt = NULL;
250 struct GMTAPI_CTRL *API = GMT->parent;
251 size_t n;
252 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
253
254 switch (opt->option) {
255 case '<': /* Input file (only one or three is accepted) */
256 Ctrl->In.active = true;
257 if (n_files >= 3) {n_errors++; continue; }
258 file[n_files] = strdup (opt->arg);
259 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(file[n_files]))) n_errors++;
260 n_files++;
261 break;
262 case '>': /* Output file (probably for -A via external interface) */
263 Ctrl->Out.active = true;
264 if (Ctrl->Out.file == NULL)
265 Ctrl->Out.file = strdup (opt->arg);
266 else
267 n_errors++;
268 break;
269
270 /* Processes program-specific parameters */
271
272 case 'A': /* Get image file name plus driver name to write via GDAL */
273 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
274 Ctrl->A.active = true;
275 if (API->external) { /* External interface only */
276 if ((n = strlen (opt->arg)) > 0) {
277 GMT_Report (API, GMT_MSG_ERROR, "Option -A: No output argument allowed\n");
278 n_errors++;
279 }
280 Ctrl->A.return_image = true;
281 Ctrl->A.way = 1; /* Building image directly, use TRPa layout, no call to GDAL */
282 }
283 else if ((n = strlen (opt->arg)) == 0) {
284 GMT_Report (API, GMT_MSG_ERROR, "Option -A: No output name provided\n");
285 n_errors++;
286 }
287 else if (gmt_M_file_is_memory (opt->arg)) {
288 Ctrl->A.file = strdup (opt->arg);
289 Ctrl->A.way = 1; /* Building image directly, use TRPa layout, no call to GDAL */
290 }
291 else if ((c = gmt_get_ext (opt->arg)) && !strcmp (c, "ppm")) { /* Want a ppm image which we can do without GDAL */
292 Ctrl->A.file = strdup (opt->arg);
293 Ctrl->A.way = 1; /* Building image directly, use TRP layout, no call to GDAL, writing a PPM file */
294 }
295 #ifdef HAVE_GDAL
296 else { /* Must give file and GDAL driver and this requires GDAL support */
297 Ctrl->A.file = strdup (opt->arg);
298 while (Ctrl->A.file[n] != '=' && n > 0) n--;
299 if (n == 0) { /* Gave no driver, see if we requested one of the standard image formats */
300 n = strlen (Ctrl->A.file) - 1;
301 while (n && Ctrl->A.file[n] != '.') n--;
302 if (n == 0) { /* Gave no image extension either... */
303 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Missing image extension or =<driver> name.\n");
304 n_errors++;
305 }
306 else { /* Check if we got a recognized extension */
307 unsigned int k, found = 0;
308 n++; /* Start of extension */
309 for (k = 0; !found && k < N_IMG_EXTENSIONS; k++) {
310 if (!strcmp (&(Ctrl->A.file[n]), gdal_ext[k]))
311 found = 1;
312 }
313 if (found == 0) {
314 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Missing the required =<driver> name.\n");
315 n_errors++;
316 }
317 else
318 GMT_Report (API, GMT_MSG_DEBUG, "grdimage: Auto-recognized GDAL driver from filename extension.\n");
319 }
320 }
321 }
322 #else
323 else
324 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Your selection requires building GMT with GDAL support.\n");
325 #endif
326 break;
327
328 case 'C': /* CPT */
329 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
330 Ctrl->C.active = true;
331 gmt_M_str_free (Ctrl->C.file);
332 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
333 gmt_cpt_interval_modifier (GMT, &(Ctrl->C.file), &(Ctrl->C.dz));
334 break;
335 #ifdef HAVE_GDAL
336 case 'D': /* Get an image via GDAL */
337 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
338 if (!API->external) /* For externals we actually still need the -D */
339 GMT_Report (API, GMT_MSG_COMPAT,
340 "Option -D is deprecated; images are detected automatically\n");
341 Ctrl->D.active = true;
342 Ctrl->D.mode = (opt->arg[0] == 'r');
343 break;
344 #endif
345 case 'E': /* Sets dpi */
346 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
347 Ctrl->E.active = true;
348 if (opt->arg[0] == 'i') /* Interpolate image to device resolution */
349 Ctrl->E.device_dpi = true;
350 else if (opt->arg[0] == '\0')
351 Ctrl->E.dpi = 100; /* Default grid dpi */
352 else
353 Ctrl->E.dpi = atoi (opt->arg);
354 break;
355 case 'G': /* -G<color>[+b|f] 1-bit fore- or background color for transparent masks (was -G[f|b]<color>) */
356 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
357 Ctrl->G.active = true;
358 if ((c = strstr (opt->arg, "+b"))) { /* Background color */
359 ind = GMT_BGD; off = GMT_FGD; k = 0; c[0] = '\0';
360 }
361 else if ((c = strstr (opt->arg, "+f"))) { /* Foreground color */
362 ind = GMT_FGD; off = GMT_BGD; k = 0; c[0] = '\0';
363 }
364 else if (gmt_M_compat_check (GMT, 5) && strchr ("bf", opt->arg[0])) { /* No modifier given so heck if first character is f or b */
365 k = 1;
366 if (opt->arg[0] == 'b') ind = GMT_BGD, off = GMT_FGD;
367 else ind = GMT_FGD, off = GMT_BGD;
368 }
369 else { /* Modern syntax where missing modifier means +f and just color is given */
370 ind = GMT_FGD; off = GMT_BGD; k = 0;
371 }
372 if (opt->arg[k] && gmt_getrgb (GMT, &opt->arg[k], Ctrl->G.rgb[ind])) {
373 gmt_rgb_syntax (GMT, 'G', " ");
374 n_errors++;
375 }
376 else
377 Ctrl->G.rgb[off][0] = -1;
378 if (c) c[0] = '+'; /* Restore */
379 break;
380 case 'I': /* Use intensity from grid or constant or auto-compute it */
381 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
382 Ctrl->I.active = true;
383 if ((c = strstr (opt->arg, "+d"))) { /* Gave +d, so derive intensities from the input grid using default settings */
384 Ctrl->I.derive = true;
385 c[0] = '\0'; /* Chop off modifier */
386 }
387 else if ((c = gmt_first_modifier (GMT, opt->arg, "amn"))) { /* Want to control how grdgradient is run */
388 unsigned int pos = 0;
389 char p[GMT_BUFSIZ] = {""};
390 Ctrl->I.derive = true;
391 while (gmt_getmodopt (GMT, 'I', c, "amn", &pos, p, &n_errors) && n_errors == 0) {
392 switch (p[0]) {
393 case 'a': gmt_M_str_free (Ctrl->I.azimuth); Ctrl->I.azimuth = strdup (&p[1]); break;
394 case 'm': gmt_M_str_free (Ctrl->I.ambient); Ctrl->I.ambient = strdup (&p[1]); break;
395 case 'n': gmt_M_str_free (Ctrl->I.method); Ctrl->I.method = strdup (&p[1]); break;
396 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
397 }
398 }
399 c[0] = '\0'; /* Chop off all modifiers so range can be determined */
400 }
401 else if (!opt->arg[0] || !strcmp (opt->arg, "+")) { /* Gave deprecated options -I or -I+ to derive intensities from the input grid using default settings */
402 Ctrl->I.derive = true;
403 if (opt->arg[0]) opt->arg[0] = '\0'; /* Remove the single + */
404 }
405 if (opt->arg[0]) { /* Gave an argument in addition to or instead of a modifier */
406 /* If given a file then it is either a intensity grid to use as is or, if I.derive is true, an alternate grid form which to derive intensities */
407 if (!gmt_access (GMT, opt->arg, R_OK)) /* Got a physical grid */
408 Ctrl->I.file = strdup (opt->arg);
409 else if (gmt_M_file_is_remote (opt->arg)) /* Got a remote grid */
410 Ctrl->I.file = strdup (opt->arg);
411 else if (opt->arg[0] && gmt_is_float (GMT, opt->arg)) { /* Looks like a constant value */
412 Ctrl->I.value = atof (opt->arg);
413 Ctrl->I.constant = true;
414 }
415 }
416 else if (!Ctrl->I.derive) {
417 GMT_Report (API, GMT_MSG_ERROR, "Option -I: Requires a valid grid file or a constant\n");
418 n_errors++;
419 }
420 if (c) c[0] = '+'; /* Restore the plus */
421 break;
422 case 'M': /* Monochrome image */
423 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
424 Ctrl->M.active = true;
425 break;
426 case 'N': /* Do not clip at map boundary */
427 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
428 Ctrl->N.active = true;
429 break;
430 case 'Q': /* PS3 colormasking */
431 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
432 Ctrl->Q.active = true;
433 if (opt->arg[0]) { /* Change input image transparency pixel color */
434 if (gmt_getrgb (GMT, opt->arg, Ctrl->Q.rgb)) { /* Change input image transparency pixel color */
435 gmt_rgb_syntax (GMT, 'Q', " ");
436 n_errors++;
437 }
438 else
439 Ctrl->Q.transp_color = true;
440 }
441 break;
442 case 'W': /* Warn if no image, usually when called from grd2kml */
443 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
444 Ctrl->W.active = true;
445 break;
446
447 default: /* Report bad options */
448 n_errors += gmt_default_error (GMT, opt->option);
449 break;
450 }
451 }
452
453 if (n_files == 3) { /* Old-style, deprecated way of plotting images via red, green, blue grids*/
454 /* We will combine these three grids into an image instead */
455 char output[GMT_VF_LEN] = {""}, cmd[GMT_LEN512] = {""};
456 GMT_Report (API, GMT_MSG_COMPAT, "Passing three grids instead of an image is deprecated. Please consider using an image instead.\n");
457 GMT_Open_VirtualFile (API, GMT_IS_IMAGE, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output);
458 for (k = 0; k < 3; k++)
459 gmt_filename_set (file[k]); /* Replace any spaces with ASCII 29 */
460 sprintf (cmd, "%s %s %s -C -N -G%s", file[0], file[1], file[2], output);
461 if (GMT_Call_Module (API, "grdmix", GMT_MODULE_CMD, cmd)) {
462 GMT_Report (API, GMT_MSG_ERROR, "Unable to combine %s/%s/%s into an image - aborting.\n", file[0], file[1], file[2]);
463 n_errors++;
464 }
465 Ctrl->In.file = strdup (output);
466 }
467 else if (n_files == 1) /* Got a single grid or image */
468 Ctrl->In.file = strdup (file[0]);
469 for (k = 0; k < 3; k++)
470 gmt_M_str_free (file[k]);
471
472 gmt_consider_current_cpt (API, &Ctrl->C.active, &(Ctrl->C.file));
473
474 if (!GMT->common.n.active && (!Ctrl->C.active || gmt_is_cpt_master (GMT, Ctrl->C.file)))
475 /* Unless user selected -n we want the default not to exceed data range on projection when we are auto-scaling a master table */
476 n_errors += gmtinit_parse_n_option (GMT, "b+c");
477
478 if (Ctrl->D.active) { /* Only OK with memory input or GDAL support */
479 if (!gmt_M_file_is_memory (Ctrl->In.file)) {
480 #ifndef HAVE_GDAL
481 GMT_Report (API, GMT_MSG_ERROR, "Option -D: Requires building GMT with GDAL support.\n");
482 n_errors++;
483 #endif
484 }
485 }
486 if (!API->external) { /* I.e, not an External interface */
487 n_errors += gmt_M_check_condition (GMT, !(n_files == 1 || n_files == 3),
488 "Must specify one (or three [deprecated]) input file(s)\n");
489 }
490 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && !Ctrl->I.constant && !Ctrl->I.file && !Ctrl->I.derive,
491 "Option -I: Must specify intensity file, value, or modifiers\n");
492 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && Ctrl->I.derive && n_files == 3,
493 "Option -I: Cannot derive intensities when r,g,b grids are given as data\n");
494 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && Ctrl->I.derive && Ctrl->D.active,
495 "Option -I: Cannot derive intensities when an image is given as data\n");
496 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && !Ctrl->E.device_dpi && Ctrl->E.dpi <= 0,
497 "Option -E: dpi must be positive\n");
498 n_errors += gmt_M_check_condition (GMT, Ctrl->G.rgb[GMT_FGD][0] < 0 && Ctrl->G.rgb[GMT_BGD][0] < 0,
499 "Option -G: Only one of fore/back-ground can be transparent for 1-bit images\n");
500 n_errors += gmt_M_check_condition (GMT, Ctrl->M.active && Ctrl->Q.active,
501 "SOption -Q: Cannot use -M when doing colormasking\n");
502 n_errors += gmt_M_check_condition (GMT, Ctrl->E.device_dpi && Ctrl->Q.active,
503 "Option -Q: Cannot use -Ei when doing colormasking\n");
504 n_errors += gmt_M_check_condition (GMT, Ctrl->A.return_image && Ctrl->Out.file == NULL,
505 "Option -A: Must provide an output filename for image\n");
506 n_errors += gmt_M_check_condition (GMT, Ctrl->A.file && Ctrl->Out.file,
507 "Option -A, -> options: Cannot provide two output files\n");
508 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->D.active,
509 "Options -C and -D options are mutually exclusive\n");
510 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
511 }
512
grdimage_clean_global_headers(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * h)513 GMT_LOCAL unsigned int grdimage_clean_global_headers (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
514 /* Problem is that many global grids are reported with misleading -R or imprecise
515 * -R or non-global -R yet nx * xinc = 360, etc. We fix these here for GMT use. */
516 unsigned int flag = 0;
517 double wasR[4], wasI[2];
518 if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* We know x is geographic longitude */
519 if (gmt_M_360_range (h->wesn[XLO], h->wesn[XHI])) /* Exact 360 range is great */
520 flag = 0;
521 else if (fabs (h->wesn[XHI] - h->wesn[XLO] - 360.0) < GMT_CONV4_LIMIT) /* Pretty close to 360, shame on manufacturer of this grid */
522 flag |= 1; /* Sloppy, we can improve this a bit */
523 else if (fabs (h->n_columns * h->inc[GMT_X] - 360.0) < GMT_CONV4_LIMIT) {
524 /* If n*xinc = 360 and previous test failed then we do not have a repeat node. These are therefore pixel grids */
525 flag |= 2; /* Global but flawed -R and registration */
526 }
527 }
528 if (gmt_M_y_is_lat (GMT, GMT_IN)) { /* We know y is geographic latitude */
529 if (gmt_M_180_range (h->wesn[YLO], h->wesn[YHI])) /* Exact 180 range is great */
530 flag |= 0;
531 else if (fabs (h->wesn[YHI] - h->wesn[YLO] - 180.0) < GMT_CONV4_LIMIT) /* Pretty close to 180, shame on manufacturer of this grid */
532 flag |= 4; /* Sloppy, we can improve this a bit */
533 else if (fabs (h->n_rows * h->inc[GMT_Y] - 180.0) < GMT_CONV4_LIMIT)
534 flag |= 8; /* Global but flawed -R and registration */
535 }
536 if (!(flag == 9 || flag == 6)) return 0; /* Only deal with mixed cases of gridline and pixel reg in lon vs lat */
537 gmt_M_memcpy (wasR, h->wesn, 4, double);
538 gmt_M_memcpy (wasI, h->inc, 2, double);
539 if ((flag & 1) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
540 h->inc[GMT_X] = 360.0 / (h->n_columns - 1); /* Get exact increment */
541 h->wesn[XHI] = h->wesn[XLO] + 360.0;
542 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global gridline-registered grid with sloppy longitude bounds.\n");
543 }
544 if ((flag & 2) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
545 h->inc[GMT_X] = 360.0 / h->n_columns; /* Get exact increment */
546 h->wesn[XHI] = h->wesn[XLO] + h->inc[GMT_X] * (h->n_columns - 1);
547 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global pixel-registered grid posing as gridline-registered with sloppy longitude bounds.\n");
548 }
549 if ((flag & 4) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
550 h->inc[GMT_Y] = 180.0 / (h->n_rows - 1); /* Get exact increment */
551 h->wesn[YLO] = -90.0;
552 h->wesn[YHI] = 90.0;
553 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global gridline-registered grid with gridline latitude bounds.\n");
554 }
555 if ((flag & 8) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
556 h->inc[GMT_Y] = 180.0 / h->n_rows; /* Get exact increment */
557 h->wesn[YLO] = -90.0 + 0.5 * h->inc[GMT_Y];
558 h->wesn[YHI] = h->wesn[YLO] + h->inc[GMT_Y] * (h->n_rows - 1);
559 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global pixel-registered grid posing as gridline-registered with sloppy latitude bounds.\n");
560 }
561 if (flag) { /* Report the before and after regions and increments */
562 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Old region: %g/%g/%g/%g with incs %g/%g\n", wasR[XLO], wasR[XHI], wasR[YLO], wasR[YHI], wasI[GMT_X], wasI[GMT_Y]);
563 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "New region: %g/%g/%g/%g with incs %g/%g\n", h->wesn[XLO], h->wesn[XHI], h->wesn[YLO], h->wesn[YHI], h->inc[GMT_X], h->inc[GMT_Y]);
564 }
565 return 1;
566 }
567
grdimage_set_proj_limits(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * r,struct GMT_GRID_HEADER * g,bool projected,unsigned int mixed)568 GMT_LOCAL void grdimage_set_proj_limits (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *r, struct GMT_GRID_HEADER *g, bool projected, unsigned int mixed) {
569 /* Sets the projected extent of the grid given the map projection.
570 * The extreme x/y coordinates are returned in r, with inc and
571 * n_columns/n_rows set accordingly. Note that some of these may change
572 * if gmt_project_init is called at a later stage */
573
574 unsigned int i, k;
575 bool all_lats = false, all_lons = false;
576 double x, y;
577 gmt_M_unused (mixed);
578
579 r->n_columns = g->n_columns; r->n_rows = g->n_rows;
580 r->registration = g->registration;
581 r->n_bands = g->n_bands;
582
583 /* By default, use entire plot region */
584
585 gmt_M_memcpy (r->wesn, GMT->current.proj.rect, 4, double);
586
587 if (GMT->current.proj.projection_GMT == GMT_GENPER && GMT->current.proj.g_width != 0.0) return;
588 if (gmt_M_is_azimuthal (GMT) && !GMT->current.proj.polar) return;
589
590 /* This fails when -R is not the entire grid region and projected is false. */
591 if (projected && gmt_M_is_geographic (GMT, GMT_IN)) {
592 all_lats = gmt_grd_is_polar (GMT, g);
593 all_lons = gmt_grd_is_global (GMT, g);
594 if (all_lons && all_lats) return; /* Whole globe */
595 }
596 /* Must search for extent along perimeter */
597
598 r->wesn[XLO] = r->wesn[YLO] = +DBL_MAX;
599 r->wesn[XHI] = r->wesn[YHI] = -DBL_MAX;
600 k = (g->registration == GMT_GRID_NODE_REG) ? 1 : 0;
601
602 for (i = 0; i < g->n_columns - k; i++) { /* South and north sides */
603 gmt_geo_to_xy (GMT, g->wesn[XLO] + i * g->inc[GMT_X], g->wesn[YLO], &x, &y);
604 r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
605 r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
606 gmt_geo_to_xy (GMT, g->wesn[XHI] - i * g->inc[GMT_X], g->wesn[YHI], &x, &y);
607 r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
608 r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
609 }
610 for (i = 0; i < g->n_rows - k; i++) { /* East and west sides */
611 gmt_geo_to_xy (GMT, g->wesn[XLO], g->wesn[YHI] - i * g->inc[GMT_Y], &x, &y);
612 r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
613 r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
614 gmt_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YLO] + i * g->inc[GMT_Y], &x, &y);
615 r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
616 r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
617 }
618 if (projected) {
619 if (all_lons) { /* Full 360, use min/max for x */
620 r->wesn[XLO] = GMT->current.proj.rect[XLO]; r->wesn[XHI] = GMT->current.proj.rect[XHI];
621 }
622 if (all_lats) { /* Full -90/+90, use min/max for y */
623 r->wesn[YLO] = GMT->current.proj.rect[YLO]; r->wesn[YHI] = GMT->current.proj.rect[YHI];
624 }
625 if (GMT->current.map.is_world && gmt_M_is_periodic (GMT)) { /* Worry about grids crossing a periodic boundary as the search above may fail */
626 if (g->wesn[XLO] < (GMT->current.proj.central_meridian+180) && g->wesn[XHI] > (GMT->current.proj.central_meridian+180))
627 r->wesn[XLO] = GMT->current.proj.rect[XLO], r->wesn[XHI] = GMT->current.proj.rect[XHI];
628 else if (g->wesn[XLO] < (GMT->current.proj.central_meridian-180) && g->wesn[XHI] > (GMT->current.proj.central_meridian-180))
629 r->wesn[XLO] = GMT->current.proj.rect[XLO], r->wesn[XHI] = GMT->current.proj.rect[XHI];
630 }
631 }
632 else if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* Extra check for non-projected longitudes that wrap */
633 double x1;
634 gmt_geo_to_xy (GMT, g->wesn[XHI]-g->inc[GMT_X], g->wesn[YHI], &x1, &y);
635 gmt_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YHI], &x, &y);
636 if (x < x1) /* Wrapped around because end of pixel is outside east; use plot width instead */
637 r->wesn[XHI] = r->wesn[XLO] + GMT->current.proj.rect[XHI];
638 }
639 }
640
grdimage_blackwhite_PS_image(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,unsigned char * image,unsigned int n_columns,unsigned int n_rows,double x_LL,double y_LL,double dx,double dy)641 GMT_LOCAL void grdimage_blackwhite_PS_image (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, unsigned char *image, unsigned int n_columns, unsigned int n_rows, double x_LL, double y_LL, double dx, double dy) {
642 /* This function takes an 8-bit grayshade image that we know is only white (255) or black (0) and converts it
643 * to a 1-bit black/white image suitable for PSL to plot using PostScripts image operator for 1-bit images.
644 * Because all projections and scalings have already taken place, this is a simple scanline operation. */
645 unsigned char *bitimage = NULL;
646 unsigned int row, col, nx8, shift, b_or_w, nx_pixels;
647 uint64_t imsize, k, k8, byte;
648 double x_side, y_side = n_rows * dy;
649
650 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Converting 8-bit image to 1-bit B/W image then plotting it\n");
651
652 nx8 = irint (ceil (n_columns / 8.0)); /* Image width must be a multiple of 8 bits, so we round up */
653 nx_pixels = nx8 * 8; /* The row length in bits after the rounding up */
654 imsize = gmt_M_get_nm (GMT, nx8, n_rows);
655 bitimage = gmt_M_memory (GMT, NULL, imsize, unsigned char); /* Memory to hold the 1-bit image */
656
657 /* Reprocess the byte image. Here there are no worries about direction of rows, cols since that was dealt with during color assignment */
658
659 for (row = 0, k = k8 = 0; row < n_rows; row++) { /* Process each scan line */
660 shift = 0; byte = 0;
661 for (col = 0; col < n_columns; col++, k++) { /* Visit each byte in the original gray shade image */
662 b_or_w = (image[k] == 255); /* Let white == 1, black == 0 */
663 byte |= b_or_w; /* Add in the current bit (0 or 1) */
664 shift++; /* Position us for the next bit */
665 if (shift == 8) { /* Did all 8, time to dump out another byte ["another one bytes the dust", he, he,...] */
666 bitimage[k8++] = (unsigned char) byte; /* Place the last 8 bits in output array... */
667 byte = shift = 0; /* ...and start over for next sequence of 8 bit nodes */
668 }
669 else /* Move the bits we have so far 1 step to the left */
670 byte <<= 1;
671 }
672 if (shift) { /* Set the remaining bits in this bit to white; this can apply to the last 1-7 nodes on each row */
673 byte |= 1;
674 shift++;
675 while (shift < 8) {
676 byte <<= 1;
677 byte |= 1;
678 shift++;
679 }
680 bitimage[k8++] = (unsigned char) byte; /* Copy final byte from this row into the image */
681 }
682 }
683
684 x_side = nx_pixels * dx; /* Since the image may be 1-7 bits wider than what we need we may enlarge it by a tiny amount */
685 PSL_plotbitimage (GMT->PSL, x_LL, y_LL, x_side, y_side, PSL_BL, bitimage, nx_pixels, n_rows, Ctrl->G.rgb[GMT_FGD], Ctrl->G.rgb[GMT_BGD]);
686 gmt_M_free (GMT, bitimage); /* Done with the B/W image array */
687 }
688
689 /* Here lies specific loop functions over grids and optional intensities etc. Because the two data sets (when intensities are involved)
690 * may very well have different padding, we must compute two node counters in some cases. When that is the case we always set these:
691 * H_s: Pointer to the header of the SOURCE grid or image
692 * H_i: Pointer to the header of the INTENSITY grid
693 * kk_s: Start index for current row of the SOURCE grid or image
694 * kk_i: Start index for current row of the intensity grid (sometimes not needed when we compute node_i directly)
695 * node_s: The current node of the SOURCE grid or image
696 * node_i: THe current node of the INTENSITY grid
697 *
698 * In the cases without intensity we simply use H_s, kk_s, and node_s for consistency.
699 */
700
grdimage_grd_gray_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)701 GMT_LOCAL void grdimage_grd_gray_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
702 /* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) no intensity */
703 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
704 uint64_t byte, kk_s, node_s;
705 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
706 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
707 gmt_M_unused (Ctrl);
708
709 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with no illumination.\n");
710 #ifdef _OPENMP
711 #pragma omp parallel for private(srow,byte,kk_s,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
712 #endif
713 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
714 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
715 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
716 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
717 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
718 (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
719 image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */
720 }
721 }
722 }
723
grdimage_grd_gray_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)724 GMT_LOCAL void grdimage_grd_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
725 /* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) with intensity */
726 int index;
727 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
728 uint64_t byte, kk_s, kk_i, node_s, node_i;
729 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
730 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
731 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
732 bool different = (Conf->int_mode && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
733
734 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with illumination.\n");
735 #ifdef _OPENMP
736 #pragma omp parallel for private(srow,byte,kk_s,kk_i,scol,node_s,node_i,index,rgb) shared(GMT,Conf,Ctrl,different,H_s,H_i,image)
737 #endif
738 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
739 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
740 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start node of this row in the data grid */
741 kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start node of same row in the intensity grid */
742 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
743 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
744 index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
745 if (index != GRDIMAGE_NAN_INDEX) { /* Add in the effect of illumination */
746 if (Conf->int_mode) { /* Intensity value comes from a co-registered grid */
747 node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
748 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
749 }
750 else /* A constant (ambient) intensity was given via -I */
751 gmt_illuminate (GMT, Ctrl->I.value, rgb);
752 }
753 image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */
754 }
755 }
756 }
757
grdimage_grd_c2s_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)758 GMT_LOCAL void grdimage_grd_c2s_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
759 /* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */
760 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
761 uint64_t byte, kk_s, node_s;
762 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
763 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
764 gmt_M_unused (Ctrl);
765
766 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
767 #ifdef _OPENMP
768 #pragma omp parallel for private(srow,byte,kk_s,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
769 #endif
770 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
771 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
772 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row in the data grid */
773 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
774 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
775 (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
776 image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
777 }
778 }
779 }
780
grdimage_grd_c2s_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)781 GMT_LOCAL void grdimage_grd_c2s_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
782 /* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */
783 int index;
784 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
785 uint64_t byte, kk_s, kk_i, node_s, node_i;
786 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
787 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
788 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
789 bool different = (Conf->int_mode && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
790
791 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
792 #ifdef _OPENMP
793 #pragma omp parallel for private(srow,byte,kk_s,kk_i,scol,node_s,node_i,index,rgb) shared(GMT,Conf,Ctrl,different,H_s,H_i,image)
794 #endif
795 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
796 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
797 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
798 kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
799 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
800 node_s = kk_s + Conf->actual_col[scol]; /* Current data grid node */
801 index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
802 if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
803 if (Conf->int_mode) { /* Intensity value comes from the grid */
804 node_i = kk_i + Conf->actual_col[scol]; /* Current intensity grid node */
805 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
806 }
807 else /* A constant (ambient) intensity was given via -I */
808 gmt_illuminate (GMT, Ctrl->I.value, rgb);
809 }
810 image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
811 }
812 }
813 }
814
grdimage_grd_color_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)815 GMT_LOCAL void grdimage_grd_color_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
816 /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
817 int k;
818 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
819 uint64_t byte, kk_s, node_s;
820 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
821 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
822 gmt_M_unused (Ctrl);
823
824 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
825 #ifdef _OPENMP
826 #pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
827 #endif
828 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
829 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
830 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
831 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
832 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
833 (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
834 for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]);
835 }
836 }
837 }
838
grdimage_grd_color_no_intensity_CM(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image,unsigned char * rgb_used)839 GMT_LOCAL void grdimage_grd_color_no_intensity_CM (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image, unsigned char *rgb_used) {
840 /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
841 unsigned char i_rgb[3];
842 int k, index;
843 int64_t srow, scol;
844 uint64_t byte, kk_s, node_s;
845 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
846 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
847 gmt_M_unused (Ctrl);
848
849 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
850 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
851 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
852 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
853 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
854 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
855 index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
856 for (k = 0; k < 3; k++) image[byte++] = i_rgb[k] = gmt_M_u255 (rgb[k]);
857 if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
858 index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */
859 rgb_used[index] = true;
860 }
861 }
862 }
863 }
864
grdimage_grd_color_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)865 GMT_LOCAL void grdimage_grd_color_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
866 /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
867 int index, k;
868 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
869 uint64_t byte, kk_s, kk_i, node_s, node_i;
870 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
871 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
872 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
873 bool different = (Conf->int_mode && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
874
875 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
876 #ifdef _OPENMP
877 #pragma omp parallel for private(srow,byte,kk_s,kk_i,k,scol,node_s,node_i,index,rgb) shared(GMT,Conf,Ctrl,different,H_s,H_i,image)
878 #endif
879 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
880 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
881 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
882 kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
883 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
884 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
885 index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
886 if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
887 if (Conf->int_mode) { /* Intensity value comes from the grid */
888 node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
889 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
890 }
891 else /* A constant (ambient) intensity was given via -I */
892 gmt_illuminate (GMT, Ctrl->I.value, rgb);
893 }
894 for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]);
895 }
896 }
897 }
898
grdimage_grd_color_with_intensity_CM(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image,unsigned char * rgb_used)899 GMT_LOCAL void grdimage_grd_color_with_intensity_CM (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image, unsigned char *rgb_used) {
900 /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
901 unsigned char i_rgb[3], n_rgb[3];
902 int index, k;
903 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
904 uint64_t byte, kk_s, kk_i, node_s, node_i;
905 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
906 struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
907 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
908 bool different = (Conf->int_mode && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens));
909
910 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
911 /* Determine NaN rgb 0-255 triple ones */
912 (void)gmt_get_rgb_from_z (GMT, Conf->P, GMT->session.d_NaN, rgb);
913 for (k = 0; k < 3; k++) n_rgb[k] = gmt_M_u255 (rgb[k]);
914
915 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
916 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
917 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
918 kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
919 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
920 node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
921 index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
922 if (index == GRDIMAGE_NAN_INDEX) { /* Nan color */
923 for (k = 0; k < 3; k++) image[byte++] = n_rgb[k];
924 }
925 else { /* Deal with illumination for non-NaN nodes */
926 if (Conf->int_mode) { /* Intensity value comes from the grid */
927 node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
928 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
929 }
930 else /* A constant (ambient) intensity was given via -I */
931 gmt_illuminate (GMT, Ctrl->I.value, rgb);
932 for (k = 0; k < 3; k++) i_rgb[k] = image[byte++] = gmt_M_u255 (rgb[k]);
933 index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */
934 rgb_used[index] = true;
935 }
936 }
937 }
938 }
939
grdimage_img_set_transparency(struct GMT_CTRL * GMT,unsigned char pix4,double * rgb)940 GMT_LOCAL void grdimage_img_set_transparency (struct GMT_CTRL *GMT, unsigned char pix4, double *rgb) {
941 /* JL: Here we assume background color is white, hence t * 1.
942 But what would it take to have a user selected background color? */
943 double o, t; /* o - opacity, t = transparency */
944 o = pix4 / 255.0; t = 1 - o;
945 rgb[0] = o * rgb[0] + t * GMT->current.map.frame.fill[GMT_Z].rgb[0];
946 rgb[1] = o * rgb[1] + t * GMT->current.map.frame.fill[GMT_Z].rgb[1];
947 rgb[2] = o * rgb[2] + t * GMT->current.map.frame.fill[GMT_Z].rgb[2];
948 }
949
grdimage_img_gray_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)950 GMT_LOCAL void grdimage_img_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
951 /* Function that fills out the image in the special case of 1) image, 2) gray, 3) with intensity */
952 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
953 uint64_t byte, kk_s, node_s, node_i;
954 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
955 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active image header */
956 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
957
958 #ifdef _OPENMP
959 #pragma omp parallel for private(srow,byte,kk_s,scol,node_s,node_i,rgb) shared(GMT,Conf,Ctrl,H_s,H_i,image)
960 #endif
961 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
962 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
963 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this image row */
964 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
965 node_s = kk_s + Conf->actual_col[scol]; /* Start of current pixel node */
966 rgb[0] = rgb[1] = rgb[2] = gmt_M_is255 (Conf->Image->data[node_s++]);
967 if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */
968 node_i = gmt_M_ijp (H_i, Conf->actual_row[srow], Conf->actual_col[scol]);
969 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb); /* Apply illumination to this color */
970 }
971 else /* A constant (ambient) intensity was given via -I */
972 gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */
973 image[byte++] = gmt_M_u255 (rgb[0]);
974 }
975 }
976 }
977
grdimage_img_gray_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)978 GMT_LOCAL void grdimage_img_gray_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
979 /* Function that fills out the image in the special case of 1) image, 2) gray, 3) no intensity */
980 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
981 uint64_t byte, kk_s, node_s;
982 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active data header */
983 gmt_M_unused (GMT);
984 gmt_M_unused (Ctrl);
985
986 #ifdef _OPENMP
987 #pragma omp parallel for private(srow,byte,kk_s,scol,node_s) shared(GMT,Conf,Ctrl,H_s,image)
988 #endif
989 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
990 byte = (uint64_t)srow * Conf->n_columns;
991 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this image row */
992 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
993 node_s = kk_s + Conf->actual_col[scol]; /* Start of current pixel node */
994 image[byte++] = Conf->Image->data[node_s];
995 }
996 }
997 }
998
grdimage_img_c2s_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)999 GMT_LOCAL void grdimage_img_c2s_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
1000 /* Function that fills out the image in the special case of 1) image, 2) color -> gray via YIQ, 3) with intensity */
1001 bool transparency = (Conf->Image->header->n_bands == 4);
1002 int64_t srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */
1003 uint64_t n_bands = Conf->Image->header->n_bands;
1004 uint64_t byte, kk_s, node_s, node_i;
1005 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
1006 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active image header */
1007 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
1008
1009 #ifdef _OPENMP
1010 #pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,node_i,rgb) shared(GMT,Conf,Ctrl,H_s,H_i,image)
1011 #endif
1012 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
1013 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
1014 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this image row */
1015 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
1016 node_s = kk_s + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */
1017 for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node_s++]);
1018 if (transparency && Conf->Image->data[node_s] < 255) /* Dealing with an image with transparency values less than 255 */
1019 grdimage_img_set_transparency (GMT, Conf->Image->data[node_s], rgb);
1020 if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */
1021 node_i = gmt_M_ijp (H_i, Conf->actual_row[srow], Conf->actual_col[scol]);
1022 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb); /* Apply illumination to this color */
1023 }
1024 else /* A constant (ambient) intensity was given via -I */
1025 gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */
1026 image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
1027 }
1028 }
1029 }
1030
grdimage_img_c2s_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)1031 GMT_LOCAL void grdimage_img_c2s_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
1032 /* Function that fills out the image in the special case of 1) image, 2) color -> gray via YIQ, 3) with intensity */
1033 bool transparency = (Conf->Image->header->n_bands == 4);
1034 int k;
1035 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
1036 uint64_t n_bands = Conf->Image->header->n_bands;
1037 uint64_t byte, kk_s, node_s;
1038 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
1039 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active data header */
1040 gmt_M_unused (Ctrl);
1041
1042 #ifdef _OPENMP
1043 #pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,rgb) shared(GMT,Conf,Ctrl,H_s,image)
1044 #endif
1045 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
1046 byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
1047 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
1048 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
1049 node_s = kk_s + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */
1050 for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node_s++]);
1051 if (transparency && Conf->Image->data[node_s] < 255) /* Dealing with an image with transparency values less than 255 */
1052 grdimage_img_set_transparency (GMT, Conf->Image->data[node_s], rgb);
1053 image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
1054 }
1055 }
1056 }
1057
grdimage_img_color_no_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)1058 GMT_LOCAL void grdimage_img_color_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
1059 /* Function that fills out the image in the special case of 1) image, 2) color, 3) with intensity */
1060 bool transparency = (Conf->Image->header->n_bands == 4);
1061 int k; /* Due to OPENMP on Windows requiring signed int loop variables */
1062 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
1063 uint64_t n_bands = Conf->Image->header->n_bands;
1064 uint64_t byte, kk_s, node_s;
1065 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
1066 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active data header */
1067 gmt_M_unused (Ctrl);
1068
1069 #ifdef _OPENMP
1070 #pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,rgb) shared(GMT,Conf,Ctrl,H_s,image)
1071 #endif
1072 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
1073 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
1074 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
1075 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
1076 node_s = kk_s + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */
1077 for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node_s++]);
1078 if (transparency && Conf->Image->data[node_s] < 255) /* Dealing with an image with transparency values less than 255 */
1079 grdimage_img_set_transparency (GMT, Conf->Image->data[node_s], rgb);
1080 for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); /* Scale up to integer 0-255 range */
1081 }
1082 }
1083 }
1084
grdimage_img_color_with_intensity(struct GMT_CTRL * GMT,struct GRDIMAGE_CTRL * Ctrl,struct GRDIMAGE_CONF * Conf,unsigned char * image)1085 GMT_LOCAL void grdimage_img_color_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
1086 /* Function that fills out the image in the special case of 1) image, 2) color, 3) with intensity */
1087 bool transparency = (Conf->Image->header->n_bands == 4);
1088 int k;
1089 int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
1090 uint64_t n_bands = Conf->Image->header->n_bands;
1091 uint64_t byte, kk_s, node_s, node_i;
1092 double rgb[4] = {0.0, 0.0, 0.0, 0.0};
1093 struct GMT_GRID_HEADER *H_s = Conf->Image->header; /* Pointer to the active image header */
1094 struct GMT_GRID_HEADER *H_i = (Conf->int_mode) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
1095
1096 #ifdef _OPENMP
1097 #pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,node_i,rgb) shared(GMT,Conf,Ctrl,H_s,H_i,image)
1098 #endif
1099 for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */
1100 byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
1101 kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this image row */
1102 for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
1103 node_s = kk_s + Conf->actual_col[scol] * n_bands; /* Start of current input pixel node */
1104 for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node_s++]);
1105 if (transparency && Conf->Image->data[node_s] < 255) /* Dealing with an image with transparency values less than 255 */
1106 grdimage_img_set_transparency (GMT, Conf->Image->data[node_s], rgb);
1107 if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */
1108 node_i = gmt_M_ijp (H_i, Conf->actual_row[srow], Conf->actual_col[scol]);
1109 gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb); /* Apply illumination to this color */
1110 }
1111 else /* A constant (ambient) intensity was given via -I */
1112 gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */
1113 for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); /* Scale up to integer 0-255 range */
1114 }
1115 }
1116 }
1117
grdimage_adjust_R_consideration(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * h)1118 GMT_LOCAL bool grdimage_adjust_R_consideration (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
1119 /* As per https://github.com/GenericMappingTools/gmt/issues/4440, when the user wants
1120 * to plot a pixel-registered global grid using a lon-lat scaling with periodic boundaries
1121 * and a central meridian that is not a multiple of the grid increment, we must actually
1122 * adjust the plot domain -R to be such a multiple. That, or require projection. */
1123
1124 double delta;
1125
1126 if (!gmt_M_is_geographic (GMT, GMT_IN)) return false; /* No geographic */
1127 if (h->registration == GMT_GRID_NODE_REG) return false; /* gridline-registration has the repeated column needed */
1128 if (gmt_M_is_nonlinear_graticule (GMT)) return true; /* Always have to project when given most projections except -JQ */
1129 if (!gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI])) return false; /* No repeating columns would be visible */
1130 delta = remainder (h->wesn[XLO] - GMT->common.R.wesn[XLO], h->inc[GMT_X]);
1131 if (gmt_M_is_zero (delta)) return false; /* No need to project if it is lining up */
1132 /* Here we need to adjust plot region */
1133 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Your grid is pixel-registered, the projection is simply longlat, and plot region is a full 360 degrees in longitude.\n");
1134 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Due to lack of repeated boundary nodes, -Rw/e must be given in multiples of the grid increment (%g)\n", h->inc[GMT_X]);
1135 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Current region (-R) setting in longitude is %g to %g\n", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
1136 GMT->common.R.wesn[XLO] += delta;
1137 GMT->common.R.wesn[XHI] += delta;
1138 if (GMT->common.R.wesn[XHI] <= 0.0) {
1139 GMT->common.R.wesn[XLO] += 360.0;
1140 GMT->common.R.wesn[XHI] += 360.0;
1141 }
1142 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Adjusted region (-R) setting in longitude is %g to %g\n", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
1143 return false;
1144 }
1145
1146 EXTERN_MSC int gmtlib_ind2rgb (struct GMT_CTRL *GMT, struct GMT_IMAGE **I_in);
1147
1148 #define bailout(code) {gmt_M_free_options (mode); return (code);}
1149 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_M_free (GMT, Conf); gmt_end_module (GMT, GMT_cpy); bailout (code);}
1150
1151 /* Two macros used to detect if an image has no meta-data for region and increments */
1152 #define img_inc_is_one(h) (h->inc[GMT_X] == 1.0 && h->inc[GMT_Y] == 1.0)
1153 #define img_region_is_dimension(h) (h->wesn[XLO] == 0.0 && h->wesn[YLO] == 0.0 && img_inc_is_one(h) && \
1154 urint (h->wesn[XHI]) == h->n_columns && urint (h->wesn[YHI]) == h->n_rows)
1155 #define img_region_is_invalid(h) (h->wesn[XLO] == 0.0 && h->wesn[YLO] == 0.0 && img_inc_is_one(h) && \
1156 (h->wesn[YHI] > 90.0 || h->wesn[XHI] > 720.0))
1157
GMT_grdimage(void * V_API,int mode,void * args)1158 EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
1159 bool done, need_to_project, normal_x, normal_y, resampled = false, gray_only = false;
1160 bool nothing_inside = false, use_intensity_grid = false, got_data_tiles = false, rgb_cube_scan;
1161 bool has_content, mem_G = false, mem_I = false, mem_D = false, got_z_grid = true;
1162 unsigned int grid_registration = GMT_GRID_NODE_REG, try, row, col, mixed = 0, pad_mode = 0;
1163 uint64_t node, k, kk, dim[GMT_DIM_SIZE] = {0, 0, 3, 0};
1164 int error = 0, ret_val = GMT_NOERROR, ftype = GMT_NOTSET;
1165
1166 char *img_ProjectionRefPROJ4 = NULL, *way[2] = {"via GDAL", "directly"}, cmd[GMT_LEN256] = {""}, data_grd[GMT_VF_LEN] = {""}, *e = NULL;
1167 unsigned char *bitimage_8 = NULL, *bitimage_24 = NULL, *rgb_used = NULL;
1168
1169 double dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0;
1170 double img_wesn[4], img_inc[2] = {1.0, 1.0}; /* Image increments & min/max for writing images or external interfaces */
1171 double *NaN_rgb = NULL, red[4] = {1.0, 0.0, 0.0, 0.0}, black[4] = {0.0, 0.0, 0.0, 0.0}, wesn[4] = {0.0, 0.0, 0.0, 0.0};
1172 double *Ix = NULL, *Iy = NULL; /* Pointers to hold on to read-only x/y image arrays */
1173
1174 struct GMT_GRID *Grid_orig = NULL, *Grid_proj = NULL;
1175 struct GMT_GRID *Intens_orig = NULL, *Intens_proj = NULL;
1176 struct GMT_GRID_HEADER_HIDDEN *HH = NULL;
1177 struct GMT_PALETTE *P = NULL;
1178 struct GRDIMAGE_CTRL *Ctrl = NULL;
1179 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT internal parameters */
1180 struct GMT_OPTION *options = NULL;
1181 struct PSL_CTRL *PSL = NULL; /* General PSL internal parameters */
1182 struct GMT_GRID_HEADER *header_work = NULL; /* Pointer to a GMT header for the image or grid */
1183 struct GMT_GRID_HEADER *header_D = NULL, *header_I = NULL, *header_G = NULL;
1184 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1185 struct GMT_IMAGE *I = NULL, *Img_proj = NULL; /* A GMT image datatype, if GDAL is used */
1186 struct GMT_IMAGE *Out = NULL; /* A GMT image datatype, if external interface is used with -A */
1187 struct GMT_GRID *G2 = NULL;
1188 struct GRDIMAGE_CONF *Conf = NULL;
1189
1190 /*----------------------- Standard module initialization and parsing ----------------------*/
1191
1192 if (API == NULL) return (GMT_NOT_A_SESSION);
1193 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
1194 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
1195
1196 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
1197
1198 /* Parse the command-line arguments */
1199
1200 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
1201 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
1202 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
1203 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
1204
1205 /*---------------------------- This is the grdimage main code ----------------------------*/
1206
1207 if ((Conf = gmt_M_memory (GMT, NULL, 1, struct GRDIMAGE_CONF)) == NULL) {
1208 GMT_Report (API, GMT_MSG_ERROR, "Failed to allocate Conf structure\n");
1209 Return (GMT_MEMORY_ERROR);
1210 }
1211
1212 gmt_grd_set_datapadding (GMT, true); /* Turn on gridpadding when reading a subset */
1213
1214 use_intensity_grid = (Ctrl->I.active && !Ctrl->I.constant); /* We want to use an intensity grid */
1215 if (Ctrl->A.file) {
1216 Ctrl->Out.file = Ctrl->A.file; Ctrl->A.file = NULL; /* Only use Out.file for writing */
1217 if ((e = gmt_get_ext (Ctrl->Out.file)) && strcmp (e, "ppm")) { /* Turn off the automatic creation of aux files by GDAL */
1218 #ifdef WIN32
1219 if (_putenv ("GDAL_PAM_ENABLED=NO"))
1220 #else
1221 if (setenv ("GDAL_PAM_ENABLED", "NO", 0))
1222 #endif
1223 GMT_Report (API, GMT_MSG_WARNING, "Unable to set GDAL_PAM_ENABLED to prevent writing of auxiliary files\n");
1224 }
1225 }
1226
1227 /* Determine if grid/image is to be projected */
1228 need_to_project = (gmt_M_is_nonlinear_graticule (GMT) || Ctrl->E.dpi > 0);
1229 if (need_to_project)
1230 GMT_Report (API, GMT_MSG_DEBUG, "Projected grid is non-orthogonal, nonlinear, or dpi was changed\n");
1231 else if (Ctrl->D.active) /* If not projecting no need for a pad */
1232 gmt_set_pad(GMT, 0);
1233 pad_mode = (need_to_project) ? GMT_GRID_NEEDS_PAD2 : 0;
1234 /* Read the illumination grid header right away so we can use its region to set that of an image (if requested) */
1235 if (use_intensity_grid) { /* Illumination grid desired */
1236 if (Ctrl->I.derive) { /* Illumination grid must be derived */
1237 /* If input grid is actually a list of data tiles then we must read that list first since this will
1238 * force all tiles to be downloaded, converted, and stitched into a single grid per -R. This must
1239 * happen _before_ we auto-derive intensities via grdgradient so that there is an input data grid */
1240
1241 if (Ctrl->I.file == NULL && gmt_file_is_tiled_list (API, Ctrl->In.file, NULL, NULL, NULL)) { /* Must read and stitch the data tiles first */
1242 if ((Grid_orig = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA|pad_mode, API->tile_wesn, Ctrl->In.file, NULL)) == NULL) { /* Get stitched grid */
1243 Return (API->error);
1244 }
1245 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Grid_orig, data_grd)) {
1246 Return (API->error);
1247 }
1248 got_data_tiles = true;
1249 }
1250 }
1251 else { /* Illumination grid present and can be read */
1252 mem_I = gmt_M_file_is_memory (Ctrl->I.file); /* Remember if the intensity grid was passed as a memory grid */
1253 GMT_Report (API, GMT_MSG_INFORMATION, "Read intensity grid header from file %s\n", Ctrl->I.file);
1254 if ((Intens_orig = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY|pad_mode, NULL, Ctrl->I.file, NULL)) == NULL) { /* Get header only */
1255 Return (API->error);
1256 }
1257 }
1258 }
1259
1260 #ifdef HAVE_GDAL
1261 if (!Ctrl->D.active && (ftype = gmt_raster_type (GMT, Ctrl->In.file, false)) == GMT_IS_IMAGE) {
1262 /* The input file is an ordinary image instead of a grid and -R may be required to use it */
1263 Ctrl->D.active = true;
1264 if (GMT->common.R.active[RSET] && !strstr (Ctrl->In.file, "earth_")) Ctrl->D.mode = true;
1265 }
1266
1267 if (!Ctrl->D.active && ftype == GMT_IS_GRID) { /* See if input could be an image of a kind that could also be a grid and we don't yet know what it is. Pass GMT_GRID_IS_IMAGE mode */
1268 if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY | GMT_GRID_IS_IMAGE|pad_mode, NULL, Ctrl->In.file, NULL)) != NULL) {
1269 gmtlib_read_grd_info (GMT, Ctrl->In.file, I->header); /* Re-read header as grid to ensure orig_datatype is set */
1270 HH = gmt_get_H_hidden (I->header); /* Get hidden structure */
1271 if (HH->orig_datatype == GMT_UCHAR || HH->orig_datatype == GMT_CHAR) Ctrl->D.active = true;
1272 /* Guess that if the image region goes from 0 to col/rol-dimensions then there is no region metadata and we want to set -Dr */
1273 if (Ctrl->D.active && GMT->common.R.active[RSET]) { /* Should we add -Dr or is it just a valid subset specification? */
1274 if (gmt_M_is_geographic (GMT, GMT_IN) && img_region_is_invalid (I->header)) Ctrl->D.mode = true;
1275 else if (img_region_is_dimension (I->header)) Ctrl->D.mode = true;
1276 }
1277 }
1278 }
1279 #endif
1280
1281 gmt_detect_oblique_region (GMT, Ctrl->In.file); /* Ensure a proper and smaller -R for oblique projections */
1282
1283 if (Ctrl->D.active) { /* Main input is a single image and not a grid */
1284 bool R_save = GMT->common.R.active[RSET];
1285 double *I_wesn = (API->got_remote_wesn) ? API->tile_wesn : NULL;
1286 if (I && !need_to_project) {
1287 gmt_M_memcpy (wesn, GMT->common.R.active[RSET] ? GMT->common.R.wesn : I->header->wesn, 4, double);
1288 need_to_project = (gmt_whole_earth (GMT, I->header->wesn, wesn)); /* Must project global images even if not needed for grids */
1289 }
1290 pad_mode = (need_to_project) ? GMT_GRID_NEEDS_PAD2 : 0;
1291 if (Ctrl->I.derive) { /* Cannot auto-derive intensities from an image */
1292 GMT_Report (API, GMT_MSG_WARNING, "Cannot derive intensities from an input image file; -I ignored\n");
1293 Ctrl->I.derive = use_intensity_grid = false;
1294 }
1295 if (use_intensity_grid && GMT->common.R.active[RSET]) {
1296 /* Make sure the region of the intensity grid and -R are in agreement within a noise threshold */
1297 double xnoise = Intens_orig->header->inc[GMT_X]*GMT_CONV4_LIMIT, ynoise = Intens_orig->header->inc[GMT_Y]*GMT_CONV4_LIMIT;
1298 if (GMT->common.R.wesn[XLO] < (Intens_orig->header->wesn[XLO]-xnoise) || GMT->common.R.wesn[XHI] > (Intens_orig->header->wesn[XHI]+xnoise) ||
1299 GMT->common.R.wesn[YLO] < (Intens_orig->header->wesn[YLO]-ynoise) || GMT->common.R.wesn[YHI] > (Intens_orig->header->wesn[YHI]+ynoise)) {
1300 GMT_Report (API, GMT_MSG_ERROR, "Requested region exceeds illumination extent\n");
1301 Return (GMT_RUNTIME_ERROR);
1302 }
1303 }
1304
1305 if (!Ctrl->D.mode && use_intensity_grid && !GMT->common.R.active[RSET]) /* Apply illumination to an image but no -R provided; use intensity domain as -R region */
1306 gmt_M_memcpy (GMT->common.R.wesn, Intens_orig->header->wesn, 4, double);
1307
1308 if (Ctrl->D.mode && GMT->common.R.active[RSET]) {
1309 /* Need to assign given -R as the image's -R, so cannot pass -R in when reading */
1310 GMT->common.R.active[RSET] = false; /* Temporarily turn off -R if given */
1311 I_wesn = NULL;
1312 }
1313 /* Read in the the entire image that is to be mapped */
1314 GMT_Report (API, GMT_MSG_INFORMATION, "Allocate memory and read image file %s\n", Ctrl->In.file);
1315 if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA | GMT_IMAGE_NO_INDEX | pad_mode, I_wesn, Ctrl->In.file, NULL)) == NULL) {
1316 Return (API->error);
1317 }
1318 GMT->common.R.active[RSET] = R_save; /* Restore -R if it was set */
1319 grid_registration = I->header->registration; /* This is presumably pixel registration since it is an image */
1320 if (grid_registration != GMT_GRID_PIXEL_REG)
1321 GMT_Report(API, GMT_MSG_INFORMATION, "Your image has gridline registration yet all images ought to be pixel registered.\n");
1322 mixed = grdimage_clean_global_headers (GMT, I->header); /* Possibly adjust flawed global -R settings in the image header */
1323 HH = gmt_get_H_hidden (I->header);
1324 if ((I->header->n_bands > 1 && strncmp (I->header->mem_layout, "BRP", 3)) || strncmp (I->header->mem_layout, "BR", 2))
1325 GMT_Report(API, GMT_MSG_INFORMATION, "The image memory layout (%s) may be of the wrong type. It should be BRPa.\n", I->header->mem_layout);
1326
1327 if (!Ctrl->D.mode && !Ctrl->I.active && !GMT->common.R.active[RSET]) /* No -R or -I were set. Use image dimensions as -R */
1328 gmt_M_memcpy (GMT->common.R.wesn, I->header->wesn, 4, double);
1329
1330 if ( (Ctrl->D.mode && GMT->common.R.active[RSET]) || (!Ctrl->D.mode && use_intensity_grid) ) { /* In either case the -R was set above or on command line */
1331 gmt_M_memcpy (I->header->wesn, GMT->common.R.wesn, 4, double);
1332 /* Get actual size of each pixel in user units */
1333 dx = gmt_M_get_inc (GMT, I->header->wesn[XLO], I->header->wesn[XHI], I->header->n_columns, I->header->registration);
1334 dy = gmt_M_get_inc (GMT, I->header->wesn[YLO], I->header->wesn[YHI], I->header->n_rows, I->header->registration);
1335 I->header->inc[GMT_X] = dx; I->header->inc[GMT_Y] = dy;
1336 HH->r_inc[GMT_X] = 1.0 / dx; /* Get inverse increments to avoid divisions later */
1337 HH->r_inc[GMT_Y] = 1.0 / dy;
1338 HH->grdtype = gmtlib_get_grdtype (GMT, GMT_IN, I->header); /* Since we now have a proper region we can determine geo/cart type */
1339 if (API->external) { /* Cannot free and update read-only image so just juggle new and old arrays */
1340 Ix = I->x; Iy = I->y; /* Keep old arrays */
1341 }
1342 else { /* Reset the grid x/y arrays */
1343 gmt_M_free (GMT, I->x); gmt_M_free (GMT, I->y);
1344 }
1345 /* Assign image coordinate arrays */
1346 I->x = gmt_grd_coord (GMT, I->header, GMT_X);
1347 I->y = gmt_grd_coord (GMT, I->header, GMT_Y);
1348 }
1349
1350 #ifdef HAVE_GDAL
1351 if (!Ctrl->A.active && gmtlib_ind2rgb(GMT, &I)) {
1352 GMT_Report (API, GMT_MSG_ERROR, "Error converting from indexed to RGB\n");
1353 Return (API->error);
1354 }
1355 #endif
1356
1357 gray_only = (I->header->n_bands == 1 && I->n_indexed_colors == 0); /* Got a grayscale image */
1358 got_z_grid = false; /* Flag that we are using a GMT_IMAGE instead of a GMT_GRID as main input */
1359
1360 if (I->header->ProjRefPROJ4 != NULL) /* We are not using this information yet, but report it under -V */
1361 GMT_Report (API, GMT_MSG_INFORMATION, "Data projection (Proj4 type)\n\t%s\n", I->header->ProjRefPROJ4);
1362
1363 header_work = I->header; /* OK, that's what what we'll use as header region to send to gmt_grd_setregion */
1364 }
1365
1366 if (!Ctrl->D.active) { /* Read the header of the input grid */
1367 mem_G = gmt_M_file_is_memory (Ctrl->In.file); /* Remember if we got the grid via a memory file or not */
1368 if (!got_data_tiles) { /* Only avoid this step if we already read a data tile bunch earlier under the I.derive == true section above */
1369 GMT_Report (API, GMT_MSG_INFORMATION, "Read header from file %s\n", Ctrl->In.file);
1370 if ((Grid_orig = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY|pad_mode, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get header only */
1371 Return (API->error);
1372 }
1373 if ((API->error = gmt_img_sanitycheck (GMT, Grid_orig->header))) { /* Used map projection on a Mercator (i.e., already a Cartesian) grid */
1374 Return (API->error);
1375 }
1376 if (!GMT->common.J.active) { /* No map projection was supplied, set default */
1377 if ((Grid_orig->header->ProjRefWKT != NULL) || (Grid_orig->header->ProjRefPROJ4 != NULL)) {
1378 static char *Jarg[2] = {"X15c", "Q15c+"}; /* The two default projections for Cartesian and Geographic */
1379 unsigned int kind = gmt_M_is_geographic (GMT, GMT_IN);
1380 gmt_parse_common_options (GMT, "J", 'J', Jarg[kind]); /* No projection specified, use linear or equidistant */
1381 GMT->common.J.active = true; /* Now parsed */
1382 }
1383 else if (GMT->current.setting.run_mode == GMT_CLASSIC) {
1384 GMT_Report (API, GMT_MSG_ERROR, "Must specify a map projection with the -J option\n");
1385 Return (GMT_PARSE_ERROR);
1386 }
1387 }
1388 }
1389 if (!Ctrl->C.active) /* Set no specific CPT so we turn -C on to use current or default CPT */
1390 Ctrl->C.active = true; /* Use default CPT (GMT->current.setting.cpt) and autostretch or under modern reuse current CPT */
1391 }
1392
1393 if (got_z_grid) header_work = Grid_orig->header; /* OK, we are in GRID mode and this was not set previously. Do it now. */
1394
1395 /* Determine what wesn to pass to gmt_map_setup */
1396
1397 if (!GMT->common.R.active[RSET] && got_z_grid) /* -R was not set so we use the grid domain */
1398 gmt_set_R_from_grd (GMT, Grid_orig->header);
1399
1400 if (Ctrl->E.dpi == 0) grdimage_adjust_R_consideration (GMT, header_work); /* Special check for global pixel-registered plots */
1401
1402 /* Initialize the projection for the selected -R -J */
1403 if (gmt_map_setup (GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
1404
1405 /* Determine the wesn to be used to read the grid file; or bail if file is outside -R */
1406
1407 if (!gmt_grd_setregion (GMT, header_work, wesn, need_to_project * GMT->common.n.interpolant))
1408 nothing_inside = true; /* Means we can return an empty plot, possibly with a basemap */
1409 else if (use_intensity_grid && !Ctrl->I.derive && !gmt_grd_setregion (GMT, Intens_orig->header, wesn, need_to_project * GMT->common.n.interpolant))
1410 nothing_inside = true;
1411
1412 if (nothing_inside) { /* Create a blank plot (or image) */
1413 GMT_Report (API, GMT_MSG_WARNING, "No grid or image inside plot domain\n");
1414 /* No grid to plot; just do an empty map and bail */
1415 /* MISSING: Action to take if -A is in effect. Need to create an empty image and return/save it */
1416 if (Ctrl->A.active) { /* Create an empty image of the right dimensions */
1417 /* Not implemented yet */
1418 GMT_Report (API, GMT_MSG_WARNING, "No image returned\n");
1419 }
1420 else { /* Plot that blank canvas */
1421 if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
1422 gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level);
1423 gmt_set_basemap_orders (GMT, GMT_BASEMAP_FRAME_AFTER, GMT_BASEMAP_GRID_AFTER, GMT_BASEMAP_ANNOT_AFTER);
1424 GMT->current.map.frame.order = GMT_BASEMAP_AFTER; /* Move to last order since only calling gmt_map_basemap once */
1425 gmt_plotcanvas (GMT); /* Fill canvas if requested */
1426 gmt_map_basemap (GMT);
1427 gmt_plane_perspective (GMT, -1, 0.0);
1428 gmt_plotend (GMT);
1429 }
1430 if (Ctrl->W.active) ret_val = GMT_IMAGE_NO_DATA; /* Flag that output image has no data - needed by grd2kml only so far */
1431 Return (ret_val);
1432 }
1433
1434 /* Here the grid/image is inside the plot domain. The same must be true of any
1435 * auto-derived intensities we may create below */
1436
1437 if (Ctrl->I.derive) { /* Auto-create intensity grid from data grid using the now determined data region */
1438 bool got_int4_grid = false;
1439 char int_grd[GMT_VF_LEN] = {""}, int4_grd[GMT_VF_LEN] = {""};
1440 double *region = (got_data_tiles) ? header_work->wesn : wesn; /* Region to pass to grdgradient */
1441 struct GMT_GRID *I_data = NULL;
1442
1443 GMT_Report (API, GMT_MSG_INFORMATION, "Derive intensity grid from data grid %s\n", (Ctrl->I.file) ? Ctrl->I.file : data_grd);
1444 /* Create a virtual file to hold the intensity grid */
1445 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, int_grd))
1446 Return (API->error);
1447 if (Ctrl->I.file) { /* Gave a file to derive from. In case it is a tiled grid we read it in here */
1448 if ((I_data = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA|pad_mode, wesn, Ctrl->I.file, NULL)) == NULL) /* Get grid data header*/
1449 Return (API->error);
1450 /* If dimensions don't match the data grid we must resample this secondary z-grid */
1451 if (Grid_orig && (I_data->header->n_columns != Grid_orig->header->n_columns || I_data->header->n_rows != Grid_orig->header->n_rows)) {
1452 char int_z_grd[GMT_VF_LEN] = {""}, *res = "gp";
1453 if (I_data->header->wesn[XLO] > region[XLO] || I_data->header->wesn[XHI] < region[XHI] || I_data->header->wesn[YLO] > region[YLO] || I_data->header->wesn[YHI] < region[YHI]) {
1454 GMT_Report (API, GMT_MSG_ERROR, "Your secondary data grid given via -I does not cover the same area as the primary grid - aborting\n");
1455 Return (GMT_GRDIO_DOMAIN_VIOLATION);
1456 }
1457 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, int_z_grd))
1458 Return (API->error);
1459 sprintf (cmd, "%s -R%.16g/%.16g/%.16g/%.16g -I%.16g/%.16g -r%c -G%s --GMT_HISTORY=readonly ",
1460 Ctrl->I.file, region[XLO], region[XHI], region[YLO], region[YHI], Grid_orig->header->inc[GMT_X], Grid_orig->header->inc[GMT_Y], res[Grid_orig->header->registration], int_z_grd);
1461 /* Call the grdsample module */
1462 GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdsample with args %s\n", cmd);
1463 if (GMT_Call_Module (API, "grdsample", GMT_MODULE_CMD, cmd))
1464 Return (API->error);
1465 /* Destroy the header we read so we can get the revised on from grdsample */
1466 if (GMT_Destroy_Data (API, &I_data) != GMT_NOERROR)
1467 Return (API->error);
1468 /* Obtain the resampled data from the virtual file */
1469 if ((I_data = GMT_Read_VirtualFile (API, int_z_grd)) == NULL)
1470 Return (API->error);
1471 }
1472 else if ((I_data = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY|pad_mode, wesn, Ctrl->I.file, I_data)) == NULL) /* Get grid data */
1473 Return (API->error);
1474 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, I_data, int4_grd))
1475 Return (API->error);
1476 got_int4_grid = true;
1477 }
1478 /* Prepare the grdgradient arguments using selected -A -N and the data region in effect. If we read in a tiled grid
1479 * then it was made with 0/360 so we must use that region in grdgradient. For non-tiles, we read the actual grid
1480 * AFTER calling gmt_mapsetup which changes the region. The tiled region remains unchanged because it is read in
1481 * all at once as it is being assembled. This fixes https://github.com/GenericMappingTools/gmt/issues/3694 */
1482 sprintf (cmd, "-G%s -A%s -N%s+a%s -R%.16g/%.16g/%.16g/%.16g --GMT_HISTORY=readonly ",
1483 int_grd, Ctrl->I.azimuth, Ctrl->I.method, Ctrl->I.ambient, region[XLO], region[XHI], region[YLO], region[YHI]);
1484 if (got_data_tiles) /* Use the virtual file we made earlier when first getting a tiled grid */
1485 strcat (cmd, data_grd);
1486 else if (got_int4_grid) /* Use the virtual file just assigned a few lines above this call */
1487 strcat (cmd, int4_grd);
1488 else { /* Default is to use the data file; we quote it in case there are spaces in the filename */
1489 gmt_filename_set (Ctrl->In.file); /* Replace any spaces with ASCII 29 */
1490 strcat (cmd, Ctrl->In.file);
1491 gmt_filename_get (Ctrl->In.file); /* Replace any ASCII 29 with spaces */
1492 }
1493 /* Call the grdgradient module */
1494 GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdgradient with args %s\n", cmd);
1495 if (GMT_Call_Module (API, "grdgradient", GMT_MODULE_CMD, cmd))
1496 Return (API->error);
1497 /* Obtain the data from the virtual file */
1498 if ((Intens_orig = GMT_Read_VirtualFile (API, int_grd)) == NULL)
1499 Return (API->error);
1500 if (got_data_tiles)
1501 GMT_Close_VirtualFile (API, data_grd);
1502 else if (got_int4_grid) /* Use the virtual file we made earlier */
1503 GMT_Close_VirtualFile (API, int4_grd);
1504 }
1505
1506 if (got_z_grid) { /* Get grid dimensions */
1507 double *region = (gmt_file_is_tiled_list (API, Ctrl->In.file, NULL, NULL, NULL)) ? API->tile_wesn : wesn; /* Make sure we get correct dimensions if tiled grids are used */
1508 Conf->n_columns = gmt_M_get_n (GMT, region[XLO], region[XHI], Grid_orig->header->inc[GMT_X], Grid_orig->header->registration);
1509 Conf->n_rows = gmt_M_get_n (GMT, region[YLO], region[YHI], Grid_orig->header->inc[GMT_Y], Grid_orig->header->registration);
1510 }
1511 else { /* Trust the image info from GDAL to make it more stable against pixel vs grid registration troubles */
1512 Conf->n_columns = I->header->n_columns;
1513 Conf->n_rows = I->header->n_rows;
1514 }
1515
1516 if (!Ctrl->A.active) { /* Initialize the PostScript plot */
1517 if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
1518 gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level);
1519 gmt_set_basemap_orders (GMT, Ctrl->N.active ? GMT_BASEMAP_FRAME_BEFORE : GMT_BASEMAP_FRAME_AFTER, GMT_BASEMAP_GRID_AFTER, GMT_BASEMAP_ANNOT_BEFORE);
1520 gmt_plotcanvas (GMT); /* Fill canvas if requested */
1521 gmt_map_basemap (GMT);
1522 if (!Ctrl->N.active) gmt_map_clip_on (GMT, GMT->session.no_rgb, 3);
1523 }
1524
1525 /* Read the grid data, possibly via a subset specified in wesn */
1526
1527 if (got_z_grid && !got_data_tiles) { /* Only skip this step if we already read a tiled grid earlier under the (I.derive == true) section */
1528 GMT_Report (API, GMT_MSG_INFORMATION, "Allocate and read data from file %s\n", Ctrl->In.file);
1529 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY|pad_mode, wesn, Ctrl->In.file, Grid_orig) == NULL) { /* Get grid data */
1530 Return (API->error);
1531 }
1532 mixed = grdimage_clean_global_headers (GMT, Grid_orig->header); /* Possibly clean up near-global headers */
1533 }
1534
1535 /* If given, get intensity grid or compute intensities (for a constant intensity) */
1536
1537 if (use_intensity_grid) { /* Illumination wanted */
1538 double *region = (gmt_file_is_tiled_list (API, Ctrl->In.file, NULL, NULL, NULL)) ? API->tile_wesn : wesn; /* Subset to pass to GMT_Read_Data if data set is tiled */
1539 GMT_Report (API, GMT_MSG_INFORMATION, "Allocates memory and read intensity file\n");
1540
1541 /* Remember, the illumination header was already read at the start of grdimage */
1542 if (!Ctrl->I.derive && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY|pad_mode, region, Ctrl->I.file, Intens_orig) == NULL) {
1543 Return (API->error); /* Failed to read the intensity grid data */
1544 }
1545 mixed = grdimage_clean_global_headers (GMT, Intens_orig->header); /* Possibly clean up near-global headers */
1546 if (got_z_grid && (Intens_orig->header->n_columns != Grid_orig->header->n_columns ||
1547 Intens_orig->header->n_rows != Grid_orig->header->n_rows)) {
1548 GMT_Report (API, GMT_MSG_ERROR, "Dimensions of intensity grid do not match that of the data grid!\n");
1549 Return (GMT_RUNTIME_ERROR);
1550 }
1551
1552 if (Ctrl->D.active && (I->header->n_columns != Intens_orig->header->n_columns || I->header->n_rows != Intens_orig->header->n_rows)) {
1553 /* Resize illumination grid to match the dimensions of the image via a call to grdsample */
1554
1555 char in_string[GMT_VF_LEN] = {""}, out_string[GMT_VF_LEN] = {""};
1556 /* Associate the intensity grid with an open virtual file - in_string will then hold the name of this input "file" */
1557 GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Intens_orig, in_string);
1558 /* Create a virtual file to hold the resampled grid - out_string then holds the name of this output "file" */
1559 GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, out_string);
1560 /* Create the command to do the resampling via the grdsample module */
1561 sprintf (cmd, "%s -G%s -I%d+/%d+ --GMT_HISTORY=readonly", in_string, out_string, (int)Conf->n_columns, (int)Conf->n_rows);
1562 GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdsample with args %s\n", cmd);
1563 if (GMT_Call_Module (GMT->parent, "grdsample", GMT_MODULE_CMD, cmd) != GMT_NOERROR) /* Do the resampling */
1564 return (API->error);
1565 /* Obtain the resmapled intensity grid from the virtual file */
1566 G2 = GMT_Read_VirtualFile (API, out_string);
1567 if (GMT_Destroy_Data (API, &Intens_orig) != GMT_NOERROR) { /* We can now delete the original intensity grid ... */
1568 Return (API->error);
1569 }
1570 Intens_orig = G2; /* ...and point to the resampled intensity grid */
1571 }
1572 }
1573
1574 if (got_z_grid && (gmt_whole_earth (GMT, Grid_orig->header->wesn, wesn) == 1))
1575 need_to_project = true; /* This can only happen if reading a remote global geographic grid and the central meridian differs from that of the projection */
1576 else if (Ctrl->D.active && (gmt_whole_earth (GMT, I->header->wesn, wesn) == 1))
1577 need_to_project = true; /* This can only happen if reading a remote global geographic image and the central meridian differs from that of the projection */
1578
1579 /* Get or calculate a color palette file */
1580
1581 has_content = (got_z_grid) ? false : true; /* Images always have content but grids may be all NaN */
1582 if (got_z_grid) { /* Got a single grid so need to convert z to color via a CPT */
1583 if (Ctrl->C.active) { /* Read a palette file */
1584 char *cpt = gmt_cpt_default (API, Ctrl->C.file, Ctrl->In.file);
1585 if ((P = gmt_get_palette (GMT, cpt, GMT_CPT_OPTIONAL, Grid_orig->header->z_min, Grid_orig->header->z_max, Ctrl->C.dz)) == NULL) {
1586 GMT_Report (API, GMT_MSG_ERROR, "Failed to read CPT %s.\n", Ctrl->C.file);
1587 gmt_free_header (API->GMT, &header_G);
1588 gmt_free_header (API->GMT, &header_I);
1589 Return (API->error); /* Well, that did not go so well... */
1590 }
1591 if (cpt) gmt_M_str_free (cpt);
1592 gray_only = (P && P->is_gray); /* Flag that we are doing a gray scale image below */
1593 Conf->P = P;
1594 if (P && P->has_pattern) GMT_Report (API, GMT_MSG_WARNING, "Patterns in CPTs will be ignored\n");
1595 }
1596 if (Ctrl->W.active) { /* Check if there are just NaNs in the grid */
1597 for (node = 0; !has_content && node < Grid_orig->header->size; node++)
1598 if (!gmt_M_is_dnan (Grid_orig->data[node])) has_content = true;
1599 }
1600 }
1601
1602 if (need_to_project) { /* Need to resample the grid or image [and intensity grid] using the specified map projection */
1603 int nx_proj = 0, ny_proj = 0;
1604 double inc[2] = {0.0, 0.0};
1605
1606 if (got_z_grid && P && P->categorical && (GMT->common.n.interpolant != BCR_NEARNEIGHBOR || GMT->common.n.antialias)) {
1607 GMT_Report (API, GMT_MSG_WARNING, "Your CPT is categorical. Enabling -nn+a to avoid interpolation across categories.\n");
1608 GMT->common.n.interpolant = BCR_NEARNEIGHBOR;
1609 GMT->common.n.antialias = false;
1610 GMT->common.n.threshold = 0.0;
1611 HH = gmt_get_H_hidden (Grid_orig->header);
1612 HH->bcr_threshold = GMT->common.n.threshold;
1613 HH->bcr_interpolant = GMT->common.n.interpolant;
1614 HH->bcr_n = 1;
1615 }
1616 if (Ctrl->E.dpi == 0) { /* Use input # of nodes as # of projected nodes */
1617 nx_proj = Conf->n_columns;
1618 ny_proj = Conf->n_rows;
1619 }
1620 if (Ctrl->D.active) { /* Must project the input image instead */
1621 GMT_Report (API, GMT_MSG_INFORMATION, "Project the input image\n");
1622 if ((Img_proj = GMT_Duplicate_Data (API, GMT_IS_IMAGE, GMT_DUPLICATE_NONE, I)) == NULL) Return (API->error); /* Just to get a header we can change */
1623 grid_registration = GMT_GRID_PIXEL_REG; /* Force pixel */
1624 grdimage_set_proj_limits (GMT, Img_proj->header, I->header, need_to_project, mixed);
1625 if (gmt_M_err_fail (GMT, gmt_project_init (GMT, Img_proj->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration),
1626 Ctrl->In.file)) Return (GMT_PROJECTION_ERROR);
1627 if (Ctrl->A.active) /* Need to set background color to white for raster images */
1628 for (k = 0; k < 3; k++) GMT->current.setting.color_patch[GMT_NAN][k] = 1.0; /* For img GDAL write use white as bg color */
1629 gmt_set_grddim (GMT, Img_proj->header); /* Recalculate projected image dimensions */
1630 if (GMT_Create_Data (API, GMT_IS_IMAGE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL, 0, 0, Img_proj) == NULL)
1631 Return (API->error); /* Failed to allocate memory for the projected image */
1632 if (gmt_img_project (GMT, I, Img_proj, false)) Return (GMT_RUNTIME_ERROR); /* Now project the image onto the projected rectangle */
1633 if (!API->external && (GMT_Destroy_Data (API, &I) != GMT_NOERROR)) { /* Free the original image now we have projected. Use Img_proj from now on */
1634 Return (API->error); /* Failed to free the image */
1635 }
1636 }
1637 else { /* Project the grid */
1638 GMT_Report (API, GMT_MSG_INFORMATION, "Project the input grid\n");
1639 if ((Grid_proj = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, Grid_orig)) == NULL)
1640 Return (API->error); /* Just to get a header we can change */
1641 /* Determine the dimensions of the projected grid */
1642 grdimage_set_proj_limits (GMT, Grid_proj->header, Grid_orig->header, need_to_project, mixed);
1643 if (grid_registration == GMT_GRID_NODE_REG) /* Force pixel if a dpi was specified, else keep as is */
1644 grid_registration = (Ctrl->E.dpi > 0) ? GMT_GRID_PIXEL_REG : Grid_orig->header->registration;
1645 if (gmt_M_err_fail (GMT, gmt_project_init (GMT, Grid_proj->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration),
1646 Ctrl->In.file)) Return (GMT_PROJECTION_ERROR);
1647 gmt_set_grddim (GMT, Grid_proj->header); /* Recalculate projected grid dimensions */
1648 if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL, 0, 0, Grid_proj) == NULL)
1649 Return (API->error); /* Failed to allocate memory for the projected grid */
1650 if (gmt_grd_project (GMT, Grid_orig, Grid_proj, false)) Return (GMT_RUNTIME_ERROR); /* Now project the grid onto the projected rectangle */
1651 if (GMT_Destroy_Data (API, &Grid_orig) != GMT_NOERROR) { /* Free the original grid now we have projected. Use Grid_proj from now on */
1652 Return (API->error); /* Failed to free the original grid */
1653 }
1654 }
1655 if (use_intensity_grid) { /* Must also project the intensity grid */
1656 GMT_Report (API, GMT_MSG_INFORMATION, "Project the intensity grid\n");
1657 if ((Intens_proj = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, Intens_orig)) == NULL) /* Just to get a header we can change */
1658 Return (API->error);
1659 if (got_z_grid) /* Use projected grid bounds as template */
1660 gmt_M_memcpy (Intens_proj->header->wesn, Grid_proj->header->wesn, 4, double);
1661 else /* Use projected image bounds as template */
1662 gmt_M_memcpy (Intens_proj->header->wesn, Img_proj->header->wesn, 4, double);
1663
1664 if (Ctrl->E.dpi == 0) { /* Use input # of nodes as # of projected nodes */
1665 nx_proj = Intens_orig->header->n_columns;
1666 ny_proj = Intens_orig->header->n_rows;
1667 }
1668 if (gmt_M_err_fail (GMT, gmt_project_init (GMT, Intens_proj->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration),
1669 Ctrl->I.file)) Return (GMT_PROJECTION_ERROR);
1670 gmt_set_grddim (GMT, Intens_proj->header); /* Recalculate projected intensity grid dimensions */
1671 if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL, 0, 0, Intens_proj) == NULL)
1672 Return (API->error); /* Failed to allocate memory for the projected intensity grid */
1673 if (gmt_grd_project (GMT, Intens_orig, Intens_proj, false)) Return (GMT_RUNTIME_ERROR); /* Now project the intensity grid onto the projected rectangle */
1674 if (GMT_Destroy_Data (API, &Intens_orig) != GMT_NOERROR) { /* Free the original intensity grid now we have projected. Use Intens_proj from now on */
1675 Return (API->error); /* Failed to free the original intensity grid */
1676 }
1677 }
1678 resampled = true; /* Yes, we did it */
1679 }
1680 else { /* Simply set the unused Grid_proj/Intens_proj pointers to point to the original unprojected Grid_orig/Intens_orig objects */
1681 struct GMT_GRID_HEADER *tmp_header = gmt_get_header (GMT);
1682 if (got_z_grid) { /* Must get a copy of the header so we can change one without affecting the other */
1683 gmt_copy_gridheader (GMT, tmp_header, Grid_orig->header);
1684 if (mem_G) { /* Need a copy of the original headers so we can restore when done */
1685 header_G = gmt_get_header (GMT);
1686 gmt_copy_gridheader (GMT, header_G, Grid_orig->header);
1687 }
1688 Grid_proj = Grid_orig;
1689 grdimage_set_proj_limits (GMT, Grid_proj->header, tmp_header, need_to_project, mixed);
1690 }
1691 if (use_intensity_grid) {
1692 if (mem_I) { /* Need a copy of the original headers so we can restore when done */
1693 header_I = gmt_get_header (GMT);
1694 gmt_copy_gridheader (GMT, header_I, Intens_orig->header);
1695 }
1696 Intens_proj = Intens_orig;
1697 }
1698 if (got_z_grid) /* Dealing with a projected grid, ensure registration is set */
1699 grid_registration = Grid_orig->header->registration;
1700 else { /* Dealing with a projected image */
1701 gmt_copy_gridheader (GMT, tmp_header, I->header);
1702 if (mem_D) { /* Need a copy of the original headers so we can restore when done */
1703 header_D = gmt_get_header (GMT);
1704 gmt_copy_gridheader (GMT, header_D, I->header);
1705 }
1706 Img_proj = I;
1707 grdimage_set_proj_limits (GMT, Img_proj->header, tmp_header, need_to_project, mixed);
1708 }
1709 gmt_free_header (API->GMT, &tmp_header);
1710 }
1711
1712 /* From here, use Grid_proj or Img_proj plus optionally Intens_proj in making the (now) Cartesian rectangular image */
1713
1714 if (got_z_grid) { /* Dealing with a projected grid, so we only have one band [z]*/
1715 Grid_proj->header->n_bands = 1;
1716 header_work = Grid_proj->header; /* Later when need to refer to the header, use this copy */
1717 }
1718 else /* Use a different reference header for the image */
1719 header_work = Img_proj->header; /* Later when need to refer to the header, use this copy */
1720
1721 /* Assign the Conf structure members */
1722 Conf->Grid = Grid_proj;
1723 Conf->Image = Img_proj;
1724 Conf->Intens = Intens_proj;
1725 Conf->n_columns = header_work->n_columns;
1726 Conf->n_rows = header_work->n_rows;
1727 Conf->int_mode = use_intensity_grid;
1728 Conf->nm = header_work->nm;
1729
1730 NaN_rgb = (P) ? P->bfn[GMT_NAN].rgb : GMT->current.setting.color_patch[GMT_NAN]; /* Determine which color represents a NaN grid node */
1731 if (got_z_grid && Ctrl->Q.active) { /* Want colormasking via the grid's NaN entries */
1732 if (gray_only) {
1733 GMT_Report (API, GMT_MSG_INFORMATION, "Your image is gray scale only but -Q requires building a 24-bit image; your image will be expanded to 24-bit.\n");
1734 gray_only = false; /* Since we cannot do 8-bit and colormasking */
1735 NaN_rgb = (Ctrl->A.active) ? black : red; /* Arbitrarily pick red for PS as the NaN color since the entire image is gray only, or black if returning an image directly */
1736 if (P) gmt_M_memcpy (P->bfn[GMT_NAN].rgb, NaN_rgb, 4, double); /* Update the NaN color entry in the CPT in memory */
1737 }
1738 if (!Ctrl->A.return_image) rgb_used = gmt_M_memory (GMT, NULL, 256*256*256, unsigned char); /* Keep track of which colors we encounter so we can find an unused unique one later */
1739 }
1740
1741 if (Ctrl->A.active) { /* We desire a raster image to be returned, not a PostScript plot */
1742 int id, k;
1743 unsigned int this_proj = GMT->current.proj.projection;
1744 char mem_layout[5] = {""}, *pch = NULL;
1745 if (Ctrl->M.active || gray_only) dim[GMT_Z] = 1; /* Only one band requested */
1746 if (!need_to_project) { /* Stick with original -R */
1747 img_wesn[XLO] = GMT->common.R.wesn[XLO]; img_wesn[XHI] = GMT->common.R.wesn[XHI];
1748 img_wesn[YHI] = GMT->common.R.wesn[YHI]; img_wesn[YLO] = GMT->common.R.wesn[YLO];
1749 }
1750 else { /* Must use the projected limits, here in meters */
1751 img_wesn[XLO] = GMT->current.proj.rect_m[XLO]; img_wesn[XHI] = GMT->current.proj.rect_m[XHI];
1752 img_wesn[YHI] = GMT->current.proj.rect_m[YHI]; img_wesn[YLO] = GMT->current.proj.rect_m[YLO];
1753 }
1754 /* Determine raster image pixel sizes in the final units */
1755 img_inc[0] = (img_wesn[XHI] - img_wesn[XLO]) / (Conf->n_columns - !grid_registration);
1756 img_inc[1] = (img_wesn[YHI] - img_wesn[YLO]) / (Conf->n_rows - !grid_registration);
1757
1758 if (grid_registration == GMT_GRID_NODE_REG) { /* Adjust domain by 1/2 pixel since SW pixel is outside the domain */
1759 img_wesn[XLO] -= 0.5 * img_inc[0]; img_wesn[XHI] += 0.5 * img_inc[0];
1760 img_wesn[YLO] -= 0.5 * img_inc[1]; img_wesn[YHI] += 0.5 * img_inc[1];
1761 }
1762 if (Ctrl->Q.active) dim[GMT_Z]++; /* Flag to remind us that we need to allocate a transparency array in GMT_Create_Data */
1763 #ifdef HAVE_GDAL
1764 if (GMT->current.gdal_read_in.O.mem_layout[0])
1765 strcpy (mem_layout, GMT->current.gdal_read_in.O.mem_layout); /* Backup current layout */
1766 else
1767 #endif
1768 gmt_strncpy (mem_layout, "TRPa", 4); /* Don't let it be empty (maybe it will screw?) */
1769 GMT_Set_Default (API, "API_IMAGE_LAYOUT", "TRPa"); /* Set grdimage's default image memory layout */
1770
1771 if ((Out = GMT_Create_Data (API, GMT_IS_IMAGE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, dim, img_wesn, img_inc, 1, 0, NULL)) == NULL) { /* Yikes, must bail */
1772 if (rgb_used) gmt_M_free (GMT, rgb_used);
1773 gmt_free_header (API->GMT, &header_G);
1774 gmt_free_header (API->GMT, &header_I);
1775 Return (API->error); /* Well, no luck with that allocation */
1776 }
1777
1778 GMT_Set_Default (API, "API_IMAGE_LAYOUT", mem_layout); /* Reset to previous memory layout */
1779
1780 HH = gmt_get_H_hidden (Out->header);
1781 if ((pch = strstr (Ctrl->Out.file, "+c")) != NULL) { /* Check if we have +c<options> for GDAL */
1782 HH->pocket = strdup (pch);
1783 pch[0] = '\0';
1784 }
1785
1786 /* See if we have valid proj info or the chosen projection has a valid PROJ4 setting */
1787 if (header_work->ProjRefWKT != NULL)
1788 Out->header->ProjRefWKT = strdup (header_work->ProjRefWKT);
1789 else if (header_work->ProjRefPROJ4 != NULL)
1790 Out->header->ProjRefPROJ4 = strdup (header_work->ProjRefPROJ4);
1791 else if (header_work->ProjRefEPSG)
1792 Out->header->ProjRefEPSG = header_work->ProjRefEPSG;
1793 else {
1794 for (k = 0, id = GMT_NOTSET; id == GMT_NOTSET && k < GMT_N_PROJ4; k++)
1795 if (GMT->current.proj.proj4[k].id == this_proj) id = k;
1796 if (id >= 0) /* Valid projection for creating world file info */
1797 img_ProjectionRefPROJ4 = gmt_export2proj4(GMT); /* Get the proj4 string */
1798 Out->header->ProjRefPROJ4 = img_ProjectionRefPROJ4;
1799 }
1800
1801 if (Ctrl->M.active || gray_only) /* Only need a byte-array to hold this image */
1802 bitimage_8 = Out->data;
1803 else /* Need 3-byte array for a 24-bit image */
1804 bitimage_24 = Out->data;
1805 }
1806 else { /* Produce a PostScript image layer */
1807 if (Ctrl->M.active || gray_only) /* Only need a byte-array to hold this image */
1808 bitimage_8 = gmt_M_memory (GMT, NULL, header_work->nm, unsigned char);
1809 else { /* Need 3-byte array for a 24-bit image, plus possibly 3 bytes for the NaN mask color */
1810 if (Ctrl->Q.active) Conf->colormask_offset = 3;
1811 bitimage_24 = gmt_M_memory (GMT, NULL, 3 * header_work->nm + Conf->colormask_offset, unsigned char);
1812 }
1813 if (Ctrl->Q.active && !(Ctrl->M.active || gray_only)) { /* Want colormasking to set a transparent pixel color */
1814 if (P) /* Use the CPT NaN color */
1815 for (k = 0; k < 3; k++) bitimage_24[k] = gmt_M_u255 (P->bfn[GMT_NAN].rgb[k]); /* Scale the NaN rgb up to 0-255 range */
1816 else if (Ctrl->Q.transp_color)
1817 for (k = 0; k < 3; k++) bitimage_24[k] = gmt_M_u255 (Ctrl->Q.rgb[k]); /* Scale the NaN rgb up to 0-255 range */
1818 /* else we default to 0 0 0 of course */
1819 }
1820 if (P && Ctrl->Q.active && !(Ctrl->M.active || gray_only)) {
1821 for (k = 0; k < 3; k++) bitimage_24[k] = gmt_M_u255 (P->bfn[GMT_NAN].rgb[k]); /* Scale the NaN rgb up to 0-255 range */
1822 }
1823 }
1824
1825 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate image pixel colors\n");
1826
1827 /* Worry about linear projections with negative scales that may reverse the orientation of the final image */
1828 normal_x = !(GMT->current.proj.projection_GMT == GMT_LINEAR && !GMT->current.proj.xyz_pos[0] && !resampled);
1829 normal_y = !(GMT->current.proj.projection_GMT == GMT_LINEAR && !GMT->current.proj.xyz_pos[1] && !resampled);
1830
1831 Conf->actual_row = gmt_M_memory (GMT, NULL, Conf->n_rows, unsigned int); /* Deal with any reversal of the y-axis due to -J */
1832 for (row = 0; row < (unsigned int)Conf->n_rows; row++) Conf->actual_row[row] = (normal_y) ? row : Conf->n_rows - row - 1;
1833 Conf->actual_col = gmt_M_memory (GMT, NULL, Conf->n_columns, unsigned int); /* Deal with any reversal of the x-axis due to -J */
1834 for (col = 0; col < (unsigned int)Conf->n_columns; col++) Conf->actual_col[col] = (normal_x) ? col : Conf->n_columns - col - 1;
1835
1836 rgb_cube_scan = (P && Ctrl->Q.active && !Ctrl->A.active); /* Need to look for unique rgb for PostScript masking */
1837
1838 /* Evaluate colors at least once (try = 0), but may do twice if -Q is active and we need to select another unique NaN color not used in the image */
1839 for (try = 0, done = false; !done && try < 2; try++) {
1840 if (got_z_grid) { /* Dealing with Grids and CPT lookup */
1841 if (Ctrl->I.active) { /* With intensity */
1842 if (gray_only) /* Grid, grayscale, with intensity, with intensity */
1843 grdimage_grd_gray_with_intensity (GMT, Ctrl, Conf, bitimage_8);
1844 else if (Ctrl->M.active) /* Grid, color converted to gray, with intensity */
1845 grdimage_grd_c2s_with_intensity (GMT, Ctrl, Conf, bitimage_8);
1846 else if (Ctrl->Q.active) /* Must deal with colormasking */
1847 grdimage_grd_color_with_intensity_CM (GMT, Ctrl, Conf, bitimage_24, rgb_used);
1848 else /* Grid, color, with intensity, no colormasking */
1849 grdimage_grd_color_with_intensity (GMT, Ctrl, Conf, bitimage_24);
1850 }
1851 else { /* No intensity */
1852 if (gray_only) /* Grid, grayscale, no intensity */
1853 grdimage_grd_gray_no_intensity (GMT, Ctrl, Conf, bitimage_8);
1854 else if (Ctrl->M.active) /* Grid, color converted to gray, with intensity */
1855 grdimage_grd_c2s_no_intensity (GMT, Ctrl, Conf, bitimage_8);
1856 else if (Ctrl->Q.active) /* Must deal with colormasking */
1857 grdimage_grd_color_no_intensity_CM (GMT, Ctrl, Conf, bitimage_24, rgb_used);
1858 else /* Grid, color, no intensity */
1859 grdimage_grd_color_no_intensity (GMT, Ctrl, Conf, bitimage_24);
1860 }
1861 }
1862 else { /* Dealing with an image, so no CPT lookup nor colormasking */
1863 if (Ctrl->I.active) { /* With intensity */
1864 if (gray_only) /* Image, grayscale, with intensity, with intensity */
1865 grdimage_img_gray_with_intensity (GMT, Ctrl, Conf, bitimage_8);
1866 else if (Ctrl->M.active) /* Image, color converted to gray, with intensity */
1867 grdimage_img_c2s_with_intensity (GMT, Ctrl, Conf, bitimage_8);
1868 else /* Grid, color, with intensity */
1869 grdimage_img_color_with_intensity (GMT, Ctrl, Conf, bitimage_24);
1870 }
1871 else { /* No intensity */
1872 if (gray_only) /* Image, grayscale, no intensity */
1873 grdimage_img_gray_no_intensity (GMT, Ctrl, Conf, bitimage_8);
1874 else if (Ctrl->M.active) /* Image, color converted to gray, with intensity */
1875 grdimage_img_c2s_no_intensity (GMT, Ctrl, Conf, bitimage_8);
1876 else /* Image, color, no intensity */
1877 grdimage_img_color_no_intensity (GMT, Ctrl, Conf, bitimage_24);
1878 }
1879 }
1880 if (rgb_cube_scan) { /* Fill in the RGB cube use */
1881 int index = 0, ks;
1882 /* Check that we found an unused r/g/b value so that colormasking will work as advertised */
1883 index = (gmt_M_u255(P->bfn[GMT_NAN].rgb[0])*256 + gmt_M_u255(P->bfn[GMT_NAN].rgb[1]))*256 + gmt_M_u255(P->bfn[GMT_NAN].rgb[2]); /* The index into the cube for the selected NaN color */
1884 if (rgb_used[index]) { /* This r/g/b already appears in the image as a non-NaN color; we must find a replacement NaN color */
1885 for (index = 0, ks = GMT_NOTSET; ks == GMT_NOTSET && index < 256*256*256; index++) /* Examine the entire cube */
1886 if (!rgb_used[index]) ks = index; /* Use this unused entry instead */
1887 if (ks == GMT_NOTSET) { /* This is really really unlikely, meaning image uses all 256^3 colors */
1888 GMT_Report (API, GMT_MSG_WARNING, "Colormasking will fail as there is no unused color that can represent transparency\n");
1889 done = true;
1890 }
1891 else { /* Pick the first unused color (i.e., ks) and let it play the role of the NaN color for transparency */
1892 bitimage_24[0] = (unsigned char)(ks >> 16);
1893 bitimage_24[1] = (unsigned char)((ks >> 8) & 255);
1894 bitimage_24[2] = (unsigned char)(ks & 255);
1895 GMT_Report (API, GMT_MSG_INFORMATION, "Transparency color reset from %s to color %d/%d/%d\n",
1896 gmt_putrgb (GMT, P->bfn[GMT_NAN].rgb), (int)bitimage_24[0], (int)bitimage_24[1], (int)bitimage_24[2]);
1897 for (k = 0; k < 3; k++) P->bfn[GMT_NAN].rgb[k] = gmt_M_is255 (bitimage_24[k]); /* Set new NaN color */
1898 }
1899 }
1900 else
1901 done = true; /* Original NaN color was never used by other nodes */
1902 }
1903 else
1904 done = true; /* Only doing the loop once here since no -Q */
1905 }
1906
1907 if (rgb_used) gmt_M_free (GMT, rgb_used); /* Done using the r/g/b cube */
1908 gmt_M_free (GMT, Conf->actual_row);
1909 gmt_M_free (GMT, Conf->actual_col);
1910
1911 if (use_intensity_grid) { /* Also done with the intensity grid */
1912 if (need_to_project || !got_z_grid) {
1913 if (GMT_Destroy_Data (API, &Intens_proj) != GMT_NOERROR)
1914 GMT_Report (API, GMT_MSG_ERROR, "Failed to free Intens_proj\n");
1915 }
1916 }
1917
1918 /* Get actual plot size of each pixel */
1919 dx = gmt_M_get_inc (GMT, header_work->wesn[XLO], header_work->wesn[XHI], Conf->n_columns, header_work->registration);
1920 dy = gmt_M_get_inc (GMT, header_work->wesn[YLO], header_work->wesn[YHI], Conf->n_rows, header_work->registration);
1921
1922 /* Set lower left position of image on map; this depends on grid or image registration */
1923
1924 x0 = header_work->wesn[XLO]; y0 = header_work->wesn[YLO];
1925 if (grid_registration == GMT_GRID_NODE_REG) { /* Grid registration, move 1/2 pixel down/left */
1926 x0 -= 0.5 * dx;
1927 y0 -= 0.5 * dy;
1928 }
1929
1930 /* Full rectangular dimension of the projected image in inches */
1931 x_side = dx * Conf->n_columns;
1932 y_side = dy * Conf->n_rows;
1933
1934 if (P && gray_only && !Ctrl->A.active) /* Determine if the gray image in fact is just black & white since the PostScript can then simplify */
1935 for (kk = 0, P->is_bw = true; P->is_bw && kk < header_work->nm; kk++)
1936 if (!(bitimage_8[kk] == 0 || bitimage_8[kk] == 255)) P->is_bw = false;
1937
1938 if (P && P->is_bw && !Ctrl->A.active) /* Can get away with a 1-bit image, but we must pack the original byte to 8 image bits, unless we are returning the image (8 or 24 bit only allowed) */
1939 grdimage_blackwhite_PS_image (GMT, Ctrl, bitimage_8, Conf->n_columns, Conf->n_rows, x0, y0, dx, dy);
1940 else if (gray_only || Ctrl->M.active) { /* Here we have a 1-layer 8 bit grayshade image */
1941 if (Ctrl->A.active) { /* Creating a raster image, not PostScript */
1942 GMT_Report (API, GMT_MSG_INFORMATION, "Writing 8-bit grayshade image %s\n", way[Ctrl->A.way]);
1943 if (GMT_Write_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Out.file, Out) != GMT_NOERROR)
1944 Return (API->error);
1945 }
1946 else { /* Lay down a PostScript 8-bit image */
1947 GMT_Report (API, GMT_MSG_INFORMATION, "Plotting 8-bit grayshade image\n");
1948 PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_8, Conf->n_columns, Conf->n_rows, (Ctrl->E.device_dpi ? -8 : 8));
1949 }
1950 }
1951 else { /* Dealing with a 24-bit color image */
1952 if (Ctrl->A.active) { /* Creating a raster image, not PostScript */
1953 if (Ctrl->Q.active) { /* Must initialize the transparency byte (alpha): 255 everywhere except at NaNs where it should be 0 */
1954 memset (Out->alpha, 255, header_work->nm);
1955 for (node = row = 0; row < (unsigned int)Conf->n_rows; row++) {
1956 kk = gmt_M_ijpgi (header_work, row, 0);
1957 for (col = 0; col < (unsigned int)Conf->n_columns; col++, node++) {
1958 if (gmt_M_is_fnan (Grid_proj->data[kk + col])) Out->alpha[node] = 0;
1959 }
1960 }
1961 }
1962 GMT_Report (API, GMT_MSG_INFORMATION, "Writing 24-bit color image %s\n", way[Ctrl->A.way]);
1963 if (GMT_Write_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Out.file, Out) != GMT_NOERROR)
1964 Return (API->error);
1965 }
1966 else { /* Lay down a PostScript 24-bit color image */
1967 GMT_Report (API, GMT_MSG_INFORMATION, "Plotting 24-bit color image\n");
1968 PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_24, (Ctrl->Q.active ? -1 : 1) *
1969 Conf->n_columns, Conf->n_rows, (Ctrl->E.device_dpi ? -24 : 24));
1970 }
1971 }
1972
1973 if (!Ctrl->A.active) { /* Finalize PostScript plot, possibly including basemap */
1974 if (!Ctrl->N.active) gmt_map_clip_off (GMT);
1975
1976 gmt_map_basemap (GMT);
1977 gmt_plane_perspective (GMT, -1, 0.0);
1978 gmt_plotend (GMT);
1979
1980 /* Here we are done producing the PostScript or raster image and can free the image arrays we used */
1981
1982 /* Free both bitimage arrays. gmt_M_free will not complain if they have not been used (NULL) */
1983 gmt_M_free (GMT, bitimage_8);
1984 gmt_M_free (GMT, bitimage_24);
1985 }
1986
1987 if (need_to_project && got_z_grid && GMT_Destroy_Data (API, &Grid_proj) != GMT_NOERROR) {
1988 GMT_Report (API, GMT_MSG_ERROR, "Failed to free Grid_proj\n");
1989 }
1990
1991 if (API->external && Ix && Iy) { /* Restore old arrays since we got a read-only image */
1992 gmt_M_free (GMT, I->x); gmt_M_free (GMT, I->y);
1993 I->x = Ix; I->y = Iy;
1994 }
1995
1996 if (got_z_grid) { /* If memory grid was passed in we must restore the header */
1997 if (mem_G && Grid_orig && header_G) {
1998 if (GMT->parent->object[0]->alloc_mode != GMT_ALLOC_EXTERNALLY) /* Because it would call gmt_M_str_free(to->ProjRefWKT) */
1999 gmt_copy_gridheader (GMT, Grid_orig->header, header_G);
2000 gmt_free_header (API->GMT, &header_G);
2001 }
2002 }
2003 if (mem_I && Intens_orig && header_I) { /* If memory grid was passed in we must restore the header */
2004 gmt_copy_gridheader (GMT, Intens_orig->header, header_I);
2005 gmt_free_header (API->GMT, &header_I);
2006 }
2007 if (mem_D && I && header_D) { /* If memory image was passed in we must restore the header */
2008 gmt_copy_gridheader (GMT, I->header, header_D);
2009 gmt_free_header (API->GMT, &header_D);
2010 }
2011
2012 if (Ctrl->D.active) { /* Free the projected image */
2013 if (GMT_Destroy_Data (API, &Img_proj) != GMT_NOERROR) {
2014 Return (API->error);
2015 }
2016 }
2017 if (Ctrl->C.active && GMT_Destroy_Data (API, &P) != GMT_NOERROR) {
2018 Return (API->error);
2019 }
2020
2021 /* May return a flag that the image/PS had no data (see -W), or just NO_ERROR */
2022 ret_val = (Ctrl->W.active && !has_content) ? GMT_IMAGE_NO_DATA : GMT_NOERROR;
2023 Return (ret_val);
2024 }
2025