-
Notifications
You must be signed in to change notification settings - Fork 346
/
grdimage.c
2394 lines (2203 loc) · 133 KB
/
grdimage.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2024 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* Author: Paul Wessel
* Date: 1-JAN-2010
* Version: 6 API
*
* Brief synopsis: grdimage will read a grid or image and plot the area
* with optional shading using the selected projection. Optionally, it
* may return the raw raster image.
*
*/
#include "gmt_dev.h"
#include "longopt/grdimage_inc.h"
#define THIS_MODULE_CLASSIC_NAME "grdimage"
#define THIS_MODULE_MODERN_NAME "grdimage"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Project and plot grids or images"
#define THIS_MODULE_KEYS "<G{+,CC(,IG(,>X},>IA,<ID"
#define THIS_MODULE_NEEDS "Jg"
#define THIS_MODULE_OPTIONS "->BJKOPRUVXYfnpt" GMT_OPT("Sc") GMT_ADD_x_OPT
/* These are images that GDAL knows how to read for us. */
#define N_IMG_EXTENSIONS 7
static char *gdal_ext[N_IMG_EXTENSIONS] = {"tiff", "tif", "gif", "png", "jpg", "jpeg", "bmp"};
#define GRDIMAGE_NAN_INDEX (GMT_NAN - 3)
/* Control structure for grdimage */
struct GRDIMAGE_CTRL {
struct GRDIMAGE_In { /* Input grid or image */
bool active;
char *file;
} In;
struct GRDIMAGE_Out { /* Output image under -A */
bool active;
char *file;
} Out;
struct GRDIMAGE_C { /* -C<cpt> or -C<color1>,<color2>[,<color3>,...][+i<dz>] */
bool active;
double dz; /* Rounding for min/max determined from data */
char *file;
char *savecpt; /* For when we want to save the automatically generated CPT */
} C;
struct GRDIMAGE_D { /* -D to read image instead of grid */
bool active;
bool mode; /* Use info of -R option to reference image */
} D;
struct GRDIMAGE_A { /* -A to write a raster file or return image to API */
bool active;
bool return_image;
unsigned int way;
char *file;
} A;
struct GRDIMAGE_E { /* -Ei|<dpi> */
bool active;
bool device_dpi;
unsigned int dpi;
} E;
struct GRDIMAGE_G { /* -G<rgb>[+b|f] */
bool active;
double rgb[2][4];
} G;
struct GRDIMAGE_I { /* -I[<intens_or_zfile>|<value>|<modifiers>] */
bool active;
bool constant;
bool derive;
double value;
char *azimuth; /* Default azimuth(s) for shading */
char *file;
char *method; /* Default scaling method */
char *ambient; /* Default ambient offset */
} I;
struct GRDIMAGE_M { /* -M */
bool active;
} M;
struct GRDIMAGE_N { /* -N */
bool active;
} N;
struct GRDIMAGE_Q { /* -Q[r/g/b][+i][+z<value>] */
bool active;
bool transp_color; /* true if a color was given */
bool z_given; /* true if a z-value was given */
bool mask_color; /* true if NaN color for image -Qcolor was set */
bool got_color; /* true if -Qcolor was set on an image */
bool invert; /* If true we turn transparency = 1 - transparency [i.e., opacity] */
double rgb[4]; /* Pixel value for transparency in images */
double value; /* If +z is used this z-value will give us the r/g/b via CPT */
} Q;
struct GRDIMAGE_T { /* -T[s][o[<pen>]] */
bool active;
bool skip;
bool outline;
struct GMT_PEN pen;
} T;
struct GRDIMAGE_W { /* -W */
bool active;
} W;
};
struct GRDIMAGE_TRANSP {
unsigned int alpha_count[256]; /* Transparency distribution for RGBA image */
unsigned int mode; /* value is transparency (mode == 0) or count (mode == 1) */
unsigned int value; /* Either chosen transparency (mode == 0) */
uint64_t n_transp; /* Number of different transparencies found */
uint64_t n_dominant; /* Number of top transparencies found */
};
struct GRDIMAGE_CONF {
/* Configuration structure for things to pass around to sub-functions */
bool invert; /* true if we should interpret transparencies as opacities */
unsigned int colormask_offset; /* Either 0 or 3 depending on -Q */
unsigned int int_mode; /* 2 if we are to apply illumination (grid) or 1 (constant) */
unsigned int *actual_row; /* Array with actual rows as a function of pseudo row */
unsigned int *actual_col; /* Array of actual columns as a function of pseudo col */
int64_t n_columns, n_rows; /* Signed dimensions for use in OpenMP loops (if implemented) */
uint64_t nm; /* Number of pixels in the image */
double dim[2]; /* Dimension of each pixel in plot units */
double orig[2]; /* Lower left image origin */
double tr_rgb[4]; /* Background color used to simulate image transparency [black] */
struct GMT_PALETTE *P; /* Pointer to the active palette [NULL if image] */
struct GMT_GRID *Grid; /* Pointer to the active grid [NULL if image] */
struct GMT_GRID *Intens; /* Pointer to the active intensity grid [NULL if no intensity] */
struct GMT_IMAGE *Image; /* Pointer to the active image [NULL if grid] */
struct GRDIMAGE_TRANSP *Transp; /* Pointer to the transparency structure [NULL if grid or RGB image] */
};
static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
struct GRDIMAGE_CTRL *C;
C = gmt_M_memory (GMT, NULL, 1, struct GRDIMAGE_CTRL);
/* Initialize values whose defaults are not 0/false/NULL */
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 */
C->I.azimuth = strdup ("-45.0"); /* Default azimuth for shading when -I+d is used */
C->I.method = strdup ("t1"); /* Default normalization for shading when -I+d is used */
C->I.ambient = strdup ("0"); /* Default ambient light for shading when -I+d is used */
C->Q.value = GMT->session.f_NaN; /* If -Q is used with a grid then NaNs are made transparent by default. */
C->Q.rgb[0] = C->Q.rgb[1] = C->Q.rgb[2] = 1.0; /* White transparent pixel is default if -Q is given without a color */
return (C);
}
static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *C) { /* Deallocate control structure */
if (!C) return;
gmt_M_str_free (C->In.file);
gmt_M_str_free (C->A.file);
gmt_M_str_free (C->C.file);
gmt_M_str_free (C->C.savecpt);
gmt_M_str_free (C->I.file);
gmt_M_str_free (C->I.azimuth);
gmt_M_str_free (C->I.ambient);
gmt_M_str_free (C->I.method);
gmt_M_str_free (C->Out.file);
gmt_M_free (GMT, C);
}
static int usage (struct GMTAPI_CTRL *API, int level) {
const char *A = " [-A<out_img>[=<driver>]]";
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
const char *extra[2] = {A, " [-A]"};
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
GMT_Usage (API, 0, "usage: %s %s %s%s [%s] [-C%s] [-D[r]] [-Ei|<dpi>] "
"[-G<rgb>[+b|f]] [-I[<intensgrid>|<value>|<modifiers>]] %s[-M] [-N] %s%s[-Q[<color>][+i][+z<value>]] "
"[%s] [-T[+o[<pen>]][+s]] [%s] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s]%s[%s]\n",
name, GMT_INGRID, GMT_J_OPT, extra[API->external], GMT_B_OPT, CPT_OPT_ARGS, API->K_OPT, API->O_OPT, API->P_OPT, GMT_Rgeo_OPT, GMT_U_OPT,
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);
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
gmt_ingrid_syntax (API, 0, "Name of a grid or image to plot");
GMT_Usage (API, -2, "Note: Grid z-values are in user units and will be "
"converted to rgb colors via the CPT. Alternatively, give a raster image. "
"If the image is plain (e.g., JPG, PNG, GIF) you must also give a corresponding -R.");
if (API->external) /* External interface */
GMT_Usage (API, -2, "If -D is used then <grid> is instead expected to be an image.");
GMT_Option (API, "J-");
GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
if (API->external) /* External interface */
GMT_Usage (API, 1, "\n-A Return a GMT raster image instead of a PostScript plot.");
else {
GMT_Usage (API, 1, "\n-A<out_img>[=<driver>]");
GMT_Usage (API, -2, "Save image in a raster format (.bmp, .gif, .jp[e]g, .png, .tif) instead of PostScript. "
"If filename does not have any of these extensions then append =<driver> to select "
"the desired image format. The 'driver' is the driver code name used by GDAL. "
"See GDAL documentation for available drivers. Note: any vector elements are lost.");
}
GMT_Option (API, "B-");
gmt_explain_cpt_input (API, 'C');
GMT_Usage (API, 1, "\n-D[r]");
if (API->external) /* External interface */
GMT_Usage (API, -2, "Input is an image instead of a grid. Append r to equate image region to -R region.");
else {
GMT_Usage (API, -2, "GMT automatically detects standard image formats. For non-standard formats you must "
"use -D to force it to be seen as an image. Append r to equate image region with -R region.");
}
GMT_Usage (API, 1, "\n-Ei|<dpi>");
GMT_Usage (API, -2, "Set dpi for the projected grid which must be constructed [100] "
"if -J implies a nonlinear graticule [Default gives same size as input grid]. "
"Alternatively, append i to do the interpolation in PostScript at device resolution (invalid with -Q).");
gmt_rgb_syntax (API->GMT, 'G', "Set transparency color for images that otherwise would result in 1-bit images. ");
GMT_Usage (API, 3, "+b Set background color.");
GMT_Usage (API, 3, "+f Set foreground color [Default].");
GMT_Usage (API, 1, "\n-I[<intensgrid>|<value>|<modifiers>]");
GMT_Usage (API, -2, "Apply directional illumination. Append name of an intensity grid, or "
"for a constant intensity (i.e., change the ambient light), just give a scalar. "
"To derive intensities from <grid> instead, append desired modifiers:");
GMT_Usage (API, 3, "+a Specify <azim> of illumination [-45].");
GMT_Usage (API, 3, "+n Set the <method> and <scale> to use [t1].");
GMT_Usage (API, 3, "+m Set <ambient> light to add [0].");
GMT_Usage (API, -2, "Alternatively, use -I+d to accept the default values (see grdgradient for more details). "
"To derive intensities from another grid than <grid>, give the alternative data grid with suitable modifiers.");
GMT_Option (API, "K");
GMT_Usage (API, 1, "\n-M Force a monochrome (gray-scale) image.");
GMT_Usage (API, 1, "\n-N Do Not clip image at the map boundary.");
GMT_Option (API, "O,P");
GMT_Usage (API, 1, "\n-Q[<color>][+i][+z<value>]");
GMT_Usage (API, -2, "Use color-masking to make grid nodes with z = NaN or white image pixels transparent. "
"Append an alternate <color> to change the transparent pixel for images [white]. "
"By default we simulate image transparency by blending pixel colors with <color> according to transparency, "
"yielding an opaque image nonetheless due to PostScript limitations. Two modifiers are available:");
GMT_Usage (API, 3, "+i Invert presumed opacities to transparencies (i.e., 1 - opacity).");
GMT_Usage (API, 3, "+z Specify grid <value> to set transparent pixel color via CPT lookup.");
GMT_Option (API, "R");
gmt_pen_syntax (API->GMT, 'T', NULL, "Image the data without interpolation by painting polygonal tiles in the form [+o[<pen>]][+s].", NULL, 0);
GMT_Usage (API, 3, "+s Skip tiles for nodes with z = NaN [Default paints all tiles].");
GMT_Usage (API, 3, "+o to draw tile outline; optionally append <pen> [Default uses no outline].");
GMT_Option (API, "U,V,X,c,f,n,p,t,x,.");
return (GMT_MODULE_USAGE);
}
/* A few non-exported library functions we need here only */
EXTERN_MSC int gmtlib_get_grdtype (struct GMT_CTRL *GMT, unsigned int direction, struct GMT_GRID_HEADER *h);
EXTERN_MSC int gmtlib_read_grd_info (struct GMT_CTRL *GMT, char *file, struct GMT_GRID_HEADER *header);
static int parse (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GMT_OPTION *options) {
/* This parses the options provided to grdimage and sets parameters in Ctrl.
* Note Ctrl has already been initialized and non-zero default values set.
* Any GMT common options will override values set previously by other commands.
* It also replaces any file names specified as input or output with the data ID
* returned when registering these sources/destinations with the API.
*/
unsigned int n_errors = 0, n_files = 0, ind, off, k;
char *c = NULL, *file[3] = {NULL, NULL, NULL}, *f = NULL;
struct GMT_OPTION *opt = NULL;
struct GMTAPI_CTRL *API = GMT->parent;
size_t n;
for (opt = options; opt; opt = opt->next) { /* Process all the options given */
switch (opt->option) {
case '<': /* Input file (only one or three is accepted) */
Ctrl->In.active = true;
if (n_files >= 3) {n_errors++; continue; }
n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(file[n_files]));
n_files++;
break;
case '>': /* Output file (probably for -A via external interface) */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Out.active);
n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_IMAGE, GMT_OUT, GMT_FILE_REMOTE, &(Ctrl->Out.file));
break;
/* Processes program-specific parameters */
case 'A': /* Get image file name plus driver name to write via GDAL */
n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
if (API->external) { /* External interface only */
if ((n = strlen (opt->arg)) > 0) {
GMT_Report (API, GMT_MSG_ERROR, "Option -A: No output argument allowed\n");
n_errors++;
}
Ctrl->A.return_image = true;
Ctrl->A.way = 1; /* Building image directly, use TRPa layout, no call to GDAL */
}
else if ((n = strlen (opt->arg)) == 0) {
GMT_Report (API, GMT_MSG_ERROR, "Option -A: No output name provided\n");
n_errors++;
}
else if (gmt_M_file_is_memory (opt->arg)) {
Ctrl->A.file = strdup (opt->arg);
Ctrl->A.way = 1; /* Building image directly, use TRPa layout, no call to GDAL */
}
else if ((c = gmt_get_ext (opt->arg)) && !strcmp (c, "ppm")) { /* Want a ppm image which we can do without GDAL */
Ctrl->A.file = strdup (opt->arg);
Ctrl->A.way = 1; /* Building image directly, use TRP layout, no call to GDAL, writing a PPM file */
}
else { /* Must give file and GDAL driver */
Ctrl->A.file = strdup (opt->arg);
while (Ctrl->A.file[n] != '=' && n > 0) n--;
if (n == 0) { /* Gave no driver, see if we requested one of the standard image formats */
n = strlen (Ctrl->A.file) - 1;
while (n && Ctrl->A.file[n] != '.') n--;
if (n == 0) { /* Gave no image extension either... */
GMT_Report (API, GMT_MSG_ERROR, "Option -A: Missing image extension or =<driver> name.\n");
n_errors++;
}
else { /* Check if we got a recognized extension */
unsigned int k, found = 0;
n++; /* Start of extension */
for (k = 0; !found && k < N_IMG_EXTENSIONS; k++) {
if (!strcmp (&(Ctrl->A.file[n]), gdal_ext[k]))
found = 1;
}
if (found == 0) {
GMT_Report (API, GMT_MSG_ERROR, "Option -A: Missing the required =<driver> name.\n");
n_errors++;
}
else
GMT_Report (API, GMT_MSG_DEBUG, "grdimage: Auto-recognized GDAL driver from filename extension.\n");
}
}
}
break;
case 'C': /* CPT */
n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
gmt_M_str_free (Ctrl->C.file);
if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
if (opt->arg[0] && (f = gmt_strrstr (Ctrl->C.file, "+s")) != NULL) { /* Filename has a +s<outname>, extract that part */
Ctrl->C.savecpt = &f[2];
f[0] = '\0'; /* Remove the +s<outname> from Ctrl->C.file */
}
gmt_cpt_interval_modifier (GMT, &(Ctrl->C.file), &(Ctrl->C.dz));
break;
case 'D': /* Get an image via GDAL */
n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
if (!API->external) /* For externals we actually still need the -D */
GMT_Report (API, GMT_MSG_COMPAT,
"Option -D is deprecated; images are detected automatically\n");
Ctrl->D.mode = (opt->arg[0] == 'r');
break;
case 'E': /* Sets dpi */
n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
if (opt->arg[0] == 'i') /* Interpolate image to device resolution */
Ctrl->E.device_dpi = true;
else if (opt->arg[0] == '\0')
Ctrl->E.dpi = 100; /* Default grid dpi */
else
Ctrl->E.dpi = atoi (opt->arg);
break;
case 'G': /* -G<color>[+b|f] 1-bit fore- or background color for transparent masks (was -G[f|b]<color>) */
n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
if ((c = strstr (opt->arg, "+b"))) { /* Background color */
ind = GMT_BGD; off = GMT_FGD; k = 0; c[0] = '\0';
}
else if ((c = strstr (opt->arg, "+f"))) { /* Foreground color */
ind = GMT_FGD; off = GMT_BGD; k = 0; c[0] = '\0';
}
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 */
k = 1;
if (opt->arg[0] == 'b') ind = GMT_BGD, off = GMT_FGD;
else ind = GMT_FGD, off = GMT_BGD;
}
else { /* Modern syntax where missing modifier means +f and just color is given */
ind = GMT_FGD; off = GMT_BGD; k = 0;
}
if (opt->arg[k] && gmt_getrgb (GMT, &opt->arg[k], Ctrl->G.rgb[ind])) {
gmt_rgb_syntax (GMT, 'G', " ");
n_errors++;
}
else
Ctrl->G.rgb[off][0] = -1;
if (c) c[0] = '+'; /* Restore */
break;
case 'I': /* Use intensity from grid or constant or auto-compute it */
n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
if ((c = strstr (opt->arg, "+d"))) { /* Gave +d, so derive intensities from the input grid using default settings */
Ctrl->I.derive = true;
c[0] = '\0'; /* Chop off modifier */
}
else if ((c = gmt_first_modifier (GMT, opt->arg, "amn"))) { /* Want to control how grdgradient is run */
unsigned int pos = 0;
char p[GMT_BUFSIZ] = {""};
Ctrl->I.derive = true;
while (gmt_getmodopt (GMT, 'I', c, "amn", &pos, p, &n_errors) && n_errors == 0) {
switch (p[0]) {
case 'a': gmt_M_str_free (Ctrl->I.azimuth); Ctrl->I.azimuth = strdup (&p[1]); break;
case 'm': gmt_M_str_free (Ctrl->I.ambient); Ctrl->I.ambient = strdup (&p[1]); break;
case 'n': gmt_M_str_free (Ctrl->I.method); Ctrl->I.method = strdup (&p[1]); break;
default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
}
}
c[0] = '\0'; /* Chop off all modifiers so range can be determined */
}
else if (!opt->arg[0] || !strcmp (opt->arg, "+")) { /* Gave deprecated options -I or -I+ to derive intensities from the input grid using default settings */
Ctrl->I.derive = true;
if (opt->arg[0]) opt->arg[0] = '\0'; /* Remove the single + */
}
if (opt->arg[0]) { /* Gave an argument in addition to or instead of a modifier */
/* 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 */
if (!gmt_access (GMT, opt->arg, R_OK)) /* Got a physical grid */
Ctrl->I.file = strdup (opt->arg);
else if (gmt_M_file_is_remote (opt->arg)) /* Got a remote grid */
Ctrl->I.file = strdup (opt->arg);
else if (opt->arg[0] && gmt_is_float (GMT, opt->arg)) { /* Looks like a constant value */
Ctrl->I.value = atof (opt->arg);
Ctrl->I.constant = true;
}
}
else if (!Ctrl->I.derive) {
GMT_Report (API, GMT_MSG_ERROR, "Option -I: Requires a valid grid file or a constant\n");
n_errors++;
}
if (c) c[0] = '+'; /* Restore the plus */
break;
case 'M': /* Monochrome image */
n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
case 'N': /* Do not clip at map boundary */
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
case 'Q': /* PS3 colormasking -Q[<color>][+i][+z<value>] */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
if (gmt_validate_modifiers (GMT, opt->arg, 'Q', "iz", GMT_MSG_ERROR)) n_errors++;
if ((c = gmt_first_modifier (GMT, opt->arg, "iz"))) { /* Got at least one modifier */
char txt[GMT_LEN256] = {""};
if (gmt_get_modifier (c, 'i', txt))
Ctrl->Q.invert = true;
if (gmt_get_modifier (c, 'z', txt)) {
if (txt[0]) {
Ctrl->Q.value = atof (txt);
Ctrl->Q.z_given = true;
}
else {
GMT_Report (API, GMT_MSG_ERROR, "Option -Q: The +z modifier requires a valid z-value\n");
n_errors++;
}
}
c[0] = '\0'; /* Chop off modifiers */
}
if (opt->arg[0]) { /* Change input image transparency pixel color from default white */
if (gmt_getrgb (GMT, opt->arg, Ctrl->Q.rgb)) { /* Change input image transparency pixel color */
gmt_rgb_syntax (GMT, 'Q', " ");
n_errors++;
}
else
Ctrl->Q.transp_color = Ctrl->Q.got_color = Ctrl->Q.mask_color = true;
}
else /* Default white transparency color */
Ctrl->Q.mask_color = true;
if (c) c[0] = '+'; /* Restore the modifier */
break;
case 'T': /* Tile plot -T[+s][+o<pen>] */
n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
if (opt->arg[0]) { /* Also got arguments */
if (strstr (opt->arg, "+s"))
Ctrl->T.skip = true;
if ((c = strstr (opt->arg, "+o"))) {
if (gmt_getpen (GMT, &c[2], &Ctrl->T.pen)) {
gmt_pen_syntax (GMT, 'T', NULL, " ", NULL, 0);
n_errors++;
}
else
Ctrl->T.outline = true;
}
}
break;
case 'W': /* Warn if no image, usually when called from grd2kml */
n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
default: /* Report bad options */
n_errors += gmt_default_option_error (GMT, opt);
break;
}
}
if (n_files == 3) { /* Old-style, deprecated way of plotting images via red, green, blue grids*/
/* We will combine these three grids into an image instead */
char output[GMT_VF_LEN] = {""}, cmd[GMT_LEN512] = {""};
GMT_Report (API, GMT_MSG_COMPAT, "Passing three grids instead of an image is deprecated. Please consider using an image instead.\n");
GMT_Open_VirtualFile (API, GMT_IS_IMAGE, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output);
for (k = 0; k < 3; k++)
gmt_filename_set (file[k]); /* Replace any spaces with ASCII 29 */
sprintf (cmd, "%s %s %s -C -N -G%s", file[0], file[1], file[2], output);
if (GMT_Call_Module (API, "grdmix", GMT_MODULE_CMD, cmd)) {
GMT_Report (API, GMT_MSG_ERROR, "Unable to combine %s/%s/%s into an image - aborting.\n", file[0], file[1], file[2]);
n_errors++;
}
Ctrl->In.file = strdup (output);
}
else if (n_files == 1) /* Got a single grid or image */
Ctrl->In.file = strdup (file[0]);
for (k = 0; k < 3; k++)
gmt_M_str_free (file[k]);
gmt_consider_current_cpt (API, &Ctrl->C.active, &(Ctrl->C.file));
if (!GMT->common.n.active && (!Ctrl->C.active || gmt_is_cpt_master (GMT, Ctrl->C.file)))
/* Unless user selected -n we want the default not to exceed data range on projection when we are auto-scaling a master table */
n_errors += gmt_parse_n_option (GMT, "c+c");
if (!API->external) { /* I.e, not an External interface */
n_errors += gmt_M_check_condition (GMT, !(n_files == 1 || n_files == 3),
"Must specify one (or three [deprecated]) input file(s)\n");
}
n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && !Ctrl->I.constant && !Ctrl->I.file && !Ctrl->I.derive,
"Option -I: Must specify intensity file, value, or modifiers\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && Ctrl->I.derive && n_files == 3,
"Option -I: Cannot derive intensities when r,g,b grids are given as data\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && Ctrl->I.derive && Ctrl->D.active,
"Option -I: Cannot derive intensities when an image is given as data\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && !Ctrl->E.device_dpi && Ctrl->E.dpi <= 0,
"Option -E: dpi must be positive\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->G.rgb[GMT_FGD][0] < 0 && Ctrl->G.rgb[GMT_BGD][0] < 0,
"Option -G: Only one of fore/back-ground can be transparent for 1-bit images\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->M.active && Ctrl->Q.active,
"Option -Q: Cannot use -M when doing colormasking\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->E.device_dpi && Ctrl->Q.active,
"Option -Q: Cannot use -Ei when doing colormasking\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->A.return_image && Ctrl->Out.file == NULL,
"Option -A: Must provide an output filename for image\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->A.file && Ctrl->Out.file,
"Option -A, -> options: Cannot provide two output files\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && (GMT->common.B.active[GMT_PRIMARY] || GMT->common.B.active[GMT_SECONDARY]),
"Option -A: Cannot draw base frame (-B) when creating just a raster image\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->Q.transp_color && Ctrl->Q.z_given,
"Option Q: Cannot both specify a r/g/b and a grid z-value\n");
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}
GMT_LOCAL unsigned int grdimage_clean_global_headers (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
/* Problem is that many global grids are reported with misleading -R or imprecise
* -R or non-global -R yet nx * xinc = 360, etc. We fix these here for GMT use. */
unsigned int flag = 0;
double wasR[4], wasI[2];
if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* We know x is geographic longitude */
if (gmt_M_360_range (h->wesn[XLO], h->wesn[XHI])) /* Exact 360 range is great */
flag = 0;
else if (fabs (h->wesn[XHI] - h->wesn[XLO] - 360.0) < GMT_CONV4_LIMIT) /* Pretty close to 360, shame on manufacturer of this grid */
flag |= 1; /* Sloppy, we can improve this a bit */
else if (fabs (h->n_columns * h->inc[GMT_X] - 360.0) < GMT_CONV4_LIMIT) {
/* If n*xinc = 360 and previous test failed then we do not have a repeat node. These are therefore pixel grids */
flag |= 2; /* Global but flawed -R and registration */
}
}
if (gmt_M_y_is_lat (GMT, GMT_IN)) { /* We know y is geographic latitude */
if (gmt_M_180_range (h->wesn[YLO], h->wesn[YHI])) /* Exact 180 range is great */
flag |= 0;
else if (fabs (h->wesn[YHI] - h->wesn[YLO] - 180.0) < GMT_CONV4_LIMIT) /* Pretty close to 180, shame on manufacturer of this grid */
flag |= 4; /* Sloppy, we can improve this a bit */
else if (fabs (h->n_rows * h->inc[GMT_Y] - 180.0) < GMT_CONV4_LIMIT)
flag |= 8; /* Global but flawed -R and registration */
}
if (!(flag == 9 || flag == 6)) return 0; /* Only deal with mixed cases of gridline and pixel reg in lon vs lat */
gmt_M_memcpy (wasR, h->wesn, 4, double);
gmt_M_memcpy (wasI, h->inc, 2, double);
if ((flag & 1) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
h->inc[GMT_X] = 360.0 / (h->n_columns - 1); /* Get exact increment */
h->wesn[XHI] = h->wesn[XLO] + 360.0;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global gridline-registered grid with sloppy longitude bounds.\n");
}
if ((flag & 2) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
h->inc[GMT_X] = 360.0 / h->n_columns; /* Get exact increment */
h->wesn[XHI] = h->wesn[XLO] + h->inc[GMT_X] * (h->n_columns - 1);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global pixel-registered grid posing as gridline-registered with sloppy longitude bounds.\n");
}
if ((flag & 4) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
h->inc[GMT_Y] = 180.0 / (h->n_rows - 1); /* Get exact increment */
h->wesn[YLO] = -90.0;
h->wesn[YHI] = 90.0;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global gridline-registered grid with gridline latitude bounds.\n");
}
if ((flag & 8) && h->registration == GMT_GRID_NODE_REG) { /* This grid needs to be treated as a pixel grid */
h->inc[GMT_Y] = 180.0 / h->n_rows; /* Get exact increment */
h->wesn[YLO] = -90.0 + 0.5 * h->inc[GMT_Y];
h->wesn[YHI] = h->wesn[YLO] + h->inc[GMT_Y] * (h->n_rows - 1);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected geographic global pixel-registered grid posing as gridline-registered with sloppy latitude bounds.\n");
}
if (flag) { /* Report the before and after regions and increments */
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]);
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]);
}
return 1;
}
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) {
/* Sets the projected extent of the grid given the map projection.
* The extreme x/y coordinates are returned in r, with inc and
* n_columns/n_rows set accordingly. Note that some of these may change
* if gmt_project_init is called at a later stage */
unsigned int i, k;
bool all_lats = false, all_lons = false;
double x, y;
gmt_M_unused (mixed);
r->n_columns = g->n_columns; r->n_rows = g->n_rows;
r->registration = g->registration;
r->n_bands = g->n_bands;
/* By default, use entire plot region */
gmt_M_memcpy (r->wesn, GMT->current.proj.rect, 4, double);
if (GMT->current.proj.projection_GMT == GMT_GENPER && GMT->current.proj.g_width != 0.0) return;
if (gmt_M_is_azimuthal (GMT) && !GMT->current.proj.polar) return;
/* This fails when -R is not the entire grid region and projected is false. */
if (projected && gmt_M_is_geographic (GMT, GMT_IN)) {
all_lats = gmt_grd_is_polar (GMT, g);
all_lons = gmt_grd_is_global (GMT, g);
if (all_lons && all_lats) return; /* Whole globe */
}
/* Must search for extent along perimeter */
r->wesn[XLO] = r->wesn[YLO] = +DBL_MAX;
r->wesn[XHI] = r->wesn[YHI] = -DBL_MAX;
k = (g->registration == GMT_GRID_NODE_REG) ? 1 : 0;
for (i = 0; i < g->n_columns - k; i++) { /* South and north sides */
gmt_geo_to_xy (GMT, g->wesn[XLO] + i * g->inc[GMT_X], g->wesn[YLO], &x, &y);
r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
gmt_geo_to_xy (GMT, g->wesn[XHI] - i * g->inc[GMT_X], g->wesn[YHI], &x, &y);
r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
}
for (i = 0; i < g->n_rows - k; i++) { /* East and west sides */
gmt_geo_to_xy (GMT, g->wesn[XLO], g->wesn[YHI] - i * g->inc[GMT_Y], &x, &y);
r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
gmt_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YLO] + i * g->inc[GMT_Y], &x, &y);
r->wesn[XLO] = MIN (r->wesn[XLO], x), r->wesn[XHI] = MAX (r->wesn[XHI], x);
r->wesn[YLO] = MIN (r->wesn[YLO], y), r->wesn[YHI] = MAX (r->wesn[YHI], y);
}
if (projected) {
if (all_lons) { /* Full 360, use min/max for x */
r->wesn[XLO] = GMT->current.proj.rect[XLO]; r->wesn[XHI] = GMT->current.proj.rect[XHI];
}
if (all_lats) { /* Full -90/+90, use min/max for y */
r->wesn[YLO] = GMT->current.proj.rect[YLO]; r->wesn[YHI] = GMT->current.proj.rect[YHI];
}
if (GMT->current.map.is_world && gmt_M_is_periodic (GMT)) { /* Worry about grids crossing a periodic boundary as the search above may fail */
if (g->wesn[XLO] < (GMT->current.proj.central_meridian+180) && g->wesn[XHI] > (GMT->current.proj.central_meridian+180))
r->wesn[XLO] = GMT->current.proj.rect[XLO], r->wesn[XHI] = GMT->current.proj.rect[XHI];
else if (g->wesn[XLO] < (GMT->current.proj.central_meridian-180) && g->wesn[XHI] > (GMT->current.proj.central_meridian-180))
r->wesn[XLO] = GMT->current.proj.rect[XLO], r->wesn[XHI] = GMT->current.proj.rect[XHI];
}
}
else if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* Extra check for non-projected longitudes that wrap */
double x1;
gmt_geo_to_xy (GMT, g->wesn[XHI]-g->inc[GMT_X], g->wesn[YHI], &x1, &y);
gmt_geo_to_xy (GMT, g->wesn[XHI], g->wesn[YHI], &x, &y);
if (x < x1) /* Wrapped around because end of pixel is outside east; use plot width instead */
r->wesn[XHI] = r->wesn[XLO] + GMT->current.proj.rect[XHI];
}
}
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) {
/* This function takes an 8-bit grayshade image that we know is only white (255) or black (0) and converts it
* to a 1-bit black/white image suitable for PSL to plot using PostScripts image operator for 1-bit images.
* Because all projections and scalings have already taken place, this is a simple scanline operation. */
unsigned char *bitimage = NULL;
unsigned int row, col, nx8, shift, b_or_w, nx_pixels;
uint64_t imsize, k, k8, byte;
double x_side, y_side = n_rows * dy;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Converting 8-bit image to 1-bit B/W image then plotting it\n");
nx8 = irint (ceil (n_columns / 8.0)); /* Image width must be a multiple of 8 bits, so we round up */
nx_pixels = nx8 * 8; /* The row length in bits after the rounding up */
imsize = gmt_M_get_nm (GMT, nx8, n_rows);
bitimage = gmt_M_memory (GMT, NULL, imsize, unsigned char); /* Memory to hold the 1-bit image */
/* Reprocess the byte image. Here there are no worries about direction of rows, cols since that was dealt with during color assignment */
for (row = 0, k = k8 = 0; row < n_rows; row++) { /* Process each scan line */
shift = 0; byte = 0;
for (col = 0; col < n_columns; col++, k++) { /* Visit each byte in the original gray shade image */
b_or_w = (image[k] == 255); /* Let white == 1, black == 0 */
byte |= b_or_w; /* Add in the current bit (0 or 1) */
shift++; /* Position us for the next bit */
if (shift == 8) { /* Did all 8, time to dump out another byte ["another one bytes the dust", he, he,...] */
bitimage[k8++] = (unsigned char) byte; /* Place the last 8 bits in output array... */
byte = shift = 0; /* ...and start over for next sequence of 8 bit nodes */
}
else /* Move the bits we have so far 1 step to the left */
byte <<= 1;
}
if (shift) { /* Set the remaining bits in this bit to white; this can apply to the last 1-7 nodes on each row */
byte |= 1;
shift++;
while (shift < 8) {
byte <<= 1;
byte |= 1;
shift++;
}
bitimage[k8++] = (unsigned char) byte; /* Copy final byte from this row into the image */
}
}
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 */
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]);
gmt_M_free (GMT, bitimage); /* Done with the B/W image array */
}
/* Here lies specific loop functions over grids and optional intensities etc. Because the two data sets (when intensities are involved)
* may very well have different padding, we must compute two node counters in some cases. When that is the case we always set these:
* H_s: Pointer to the header of the SOURCE grid or image
* H_i: Pointer to the header of the INTENSITY grid
* kk_s: Start index for current row of the SOURCE grid or image
* kk_i: Start index for current row of the intensity grid (sometimes not needed when we compute node_i directly)
* node_s: The current node of the SOURCE grid or image
* node_i: THe current node of the INTENSITY grid
*
* In the cases without intensity we simply use H_s, kk_s, and node_s for consistency.
*/
GMT_LOCAL void grdimage_grd_gray_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) no intensity */
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, node_s;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
gmt_M_unused (Ctrl);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with no illumination.\n");
#ifdef _OPENMP
#pragma omp parallel for private(srow,byte,kk_s,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
(void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */
}
}
}
GMT_LOCAL void grdimage_grd_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) with intensity */
int index;
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, kk_i, node_s, node_i;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
struct GMT_GRID_HEADER *H_i = (Conf->int_mode == 2) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
bool different = (Conf->int_mode == 2 && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with illumination.\n");
#ifdef _OPENMP
#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)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start node of this row in the data grid */
kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start node of same row in the intensity grid */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
if (index != GRDIMAGE_NAN_INDEX) { /* Add in the effect of illumination */
if (Conf->int_mode == 2) { /* Intensity value comes from a co-registered grid */
node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
}
else if (Conf->int_mode == 1) /* A constant (ambient) intensity was given via -I */
gmt_illuminate (GMT, Ctrl->I.value, rgb);
}
image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */
}
}
}
GMT_LOCAL void grdimage_grd_c2s_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, node_s;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
gmt_M_unused (Ctrl);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
#ifdef _OPENMP
#pragma omp parallel for private(srow,byte,kk_s,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row in the data grid */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
(void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
}
}
}
GMT_LOCAL void grdimage_grd_c2s_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */
int index;
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, kk_i, node_s, node_i;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
struct GMT_GRID_HEADER *H_i = (Conf->int_mode == 2) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
bool different = (Conf->int_mode == 2 && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
#ifdef _OPENMP
#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)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)srow * Conf->n_columns; /* Start of output gray image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current data grid node */
index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
if (Conf->int_mode == 2) { /* Intensity value comes from the grid */
node_i = kk_i + Conf->actual_col[scol]; /* Current intensity grid node */
gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
}
else if (Conf->int_mode == 1) /* A constant (ambient) intensity was given via -I */
gmt_illuminate (GMT, Ctrl->I.value, rgb);
}
image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb));
}
}
}
GMT_LOCAL void grdimage_grd_color_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
int k;
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, node_s;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
gmt_M_unused (Ctrl);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
#ifdef _OPENMP
#pragma omp parallel for private(srow,byte,kk_s,k,scol,node_s,rgb) shared(GMT,Conf,H_s,image)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
(void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]);
}
}
}
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) {
/* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
unsigned char i_rgb[3];
int k, index;
int64_t srow, scol;
uint64_t byte, kk_s, node_s;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
gmt_M_unused (Ctrl);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this row */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
for (k = 0; k < 3; k++) image[byte++] = i_rgb[k] = gmt_M_u255 (rgb[k]);
if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */
rgb_used[index] = true;
}
}
}
}
GMT_LOCAL void grdimage_grd_color_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
int index, k;
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, kk_i, node_s, node_i;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
struct GMT_GRID_HEADER *H_i = (Conf->int_mode == 2) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
bool different = (Conf->int_mode == 2 && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens)); /* True of they differ in padding */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
#ifdef _OPENMP
#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)
#endif
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */
if (Conf->int_mode == 2) { /* Intensity value comes from the grid */
node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
}
else if (Conf->int_mode == 1) /* A constant (ambient) intensity was given via -I */
gmt_illuminate (GMT, Ctrl->I.value, rgb);
}
for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]);
}
}
}
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) {
/* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */
unsigned char i_rgb[3], n_rgb[3];
int index, k;
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */
uint64_t byte, kk_s, kk_i, node_s, node_i;
double rgb[4] = {0.0, 0.0, 0.0, 0.0};
struct GMT_GRID_HEADER *H_s = Conf->Grid->header; /* Pointer to the active data header */
struct GMT_GRID_HEADER *H_i = (Conf->int_mode == 2) ? Conf->Intens->header : NULL; /* Pointer to the active intensity header */
bool different = (Conf->int_mode == 2 && !gmt_M_grd_same_region (GMT,Conf->Grid,Conf->Intens));
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n");
/* Determine NaN rgb 0-255 triple ones */
(void)gmt_get_rgb_from_z (GMT, Conf->P, GMT->session.d_NaN, rgb);
for (k = 0; k < 3; k++) n_rgb[k] = gmt_M_u255 (rgb[k]);
for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */
byte = (uint64_t)Conf->colormask_offset + 3 * srow * Conf->n_columns; /* Start of output color image row */
kk_s = gmt_M_ijpgi (H_s, Conf->actual_row[srow], 0); /* Start pixel of this data row */
kk_i = (different) ? gmt_M_ijpgi (H_i, Conf->actual_row[srow], 0) : kk_s; /* Start pixel of this row in the intensity grid */
for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */
node_s = kk_s + Conf->actual_col[scol]; /* Current grid node */
index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node_s], rgb);
if (index == GRDIMAGE_NAN_INDEX) { /* Nan color */
for (k = 0; k < 3; k++) image[byte++] = n_rgb[k];
}
else { /* Deal with illumination for non-NaN nodes */
if (Conf->int_mode == 2) { /* Intensity value comes from the grid */
node_i = kk_i + Conf->actual_col[scol]; /* Current intensity node */
gmt_illuminate (GMT, Conf->Intens->data[node_i], rgb);
}
else if (Conf->int_mode == 1) /* A constant (ambient) intensity was given via -I */
gmt_illuminate (GMT, Ctrl->I.value, rgb);
for (k = 0; k < 3; k++) i_rgb[k] = image[byte++] = gmt_M_u255 (rgb[k]);
index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */
rgb_used[index] = true;
}
}
}
}
GMT_LOCAL void grdimage_img_set_transparency (struct GMT_CTRL *GMT, struct GRDIMAGE_CONF *Conf, unsigned char pix4, double *rgb) {
/* JL: Here we assume Conf->tr_rgb background color is white, hence t * 1.
If user selected background color with -Q<rgb> that that is what Conf->tr_rgb contains */
/* PW Dec 2023: There are several situations we must handle:
* 1. There is no transparency, in which case this function is not called.
* 2. There are only two levels of transparency 0 or 255. In that case
* we let t = 1 and effectively select the -Q<rgb> color for transparency
* By default that is white.
* 3. If more than 2 transparencies we expect we have full variable transparency
* and most of the time we will get a blended color output.
*/
double o, t; /* o - opacity, t = transparency */
gmt_M_unused (GMT);
o = pix4 / 255.0; t = 1 - o;
if (Conf->Transp->n_transp == 2)
gmt_M_rgb_only_copy (rgb, Conf->tr_rgb);
else { /* Blend */
if (Conf->invert) t = o, o = 1 - t;
rgb[0] = o * rgb[0] + t * Conf->tr_rgb[0];
rgb[1] = o * rgb[1] + t * Conf->tr_rgb[1];
rgb[2] = o * rgb[2] + t * Conf->tr_rgb[2];
}
rgb[3] = 0;
}
GMT_LOCAL void grdimage_img_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) {
/* Function that fills out the image in the special case of 1) image, 2) gray, 3) with intensity */
int64_t srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */