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 * Brief synopsis: read a grid file and find the values which divide its range
19 * into n_cell number of quantiles.
20 *
21 * Author: W.H.F. Smith
22 * Date: 31 May 1990
23 * Version: 6 API
24 */
25
26 #include "gmt_dev.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "grdhisteq"
29 #define THIS_MODULE_MODERN_NAME "grdhisteq"
30 #define THIS_MODULE_LIB "core"
31 #define THIS_MODULE_PURPOSE "Perform histogram equalization for a grid"
32 #define THIS_MODULE_KEYS "<G{,GG},DD)"
33 #define THIS_MODULE_NEEDS ""
34 #define THIS_MODULE_OPTIONS "-RVh"
35
36 struct GRDHISTEQ_CTRL {
37 struct GRDHISTEQ_In {
38 bool active;
39 char *file;
40 } In;
41 struct GRDHISTEQ_C { /* -C<n_cells>*/
42 bool active;
43 unsigned int value;
44 } C;
45 struct GRDHISTEQ_D { /* -D[<file>] */
46 bool active;
47 char *file;
48 } D;
49 struct GRDHISTEQ_G { /* -G<file> */
50 bool active;
51 char *file;
52 } G;
53 struct GRDHISTEQ_N { /* -N[<norm>] */
54 bool active;
55 double norm;
56 } N;
57 struct GRDHISTEQ_Q { /* -Q */
58 bool active;
59 } Q;
60 };
61
62 struct INDEXED_DATA {
63 gmt_grdfloat x;
64 uint64_t i;
65 };
66
67 struct GRDHISTEQ_CELL {
68 unsigned int row;
69 gmt_grdfloat low;
70 gmt_grdfloat high;
71 };
72
73 EXTERN_MSC int gmtlib_compare_observation (const void *a, const void *b);
74
New_Ctrl(struct GMT_CTRL * GMT)75 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
76 struct GRDHISTEQ_CTRL *C = NULL;
77
78 C = gmt_M_memory (GMT, NULL, 1, struct GRDHISTEQ_CTRL);
79
80 /* Initialize values whose defaults are not 0/false/NULL */
81 C->C.value = 16;
82 return (C);
83 }
84
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDHISTEQ_CTRL * C)85 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDHISTEQ_CTRL *C) { /* Deallocate control structure */
86 if (!C) return;
87 gmt_M_str_free (C->In.file);
88 gmt_M_str_free (C->D.file);
89 gmt_M_str_free (C->G.file);
90 gmt_M_free (GMT, C);
91 }
92
usage(struct GMTAPI_CTRL * API,int level)93 static int usage (struct GMTAPI_CTRL *API, int level) {
94 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
95 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
96 GMT_Usage (API, 0, "usage: %s %s [-C<n_cells>] [-D[<table>]] [-G%s] [-N[<norm>]] [-Q] "
97 "[%s] [%s] [%s] [%s]\n", name, GMT_INGRID, GMT_OUTGRID, GMT_Rgeo_OPT, GMT_V_OPT, GMT_ho_OPT, GMT_PAR_OPT);
98
99 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
100
101 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
102 gmt_ingrid_syntax (API, 0, "Name of input grid");
103 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
104 GMT_Usage (API, 1, "\n-C<n_cells>");
105 GMT_Usage (API, -2, "Set how many cells (divisions) of data range to make [16].");
106 GMT_Usage (API, 1, "\n-D[<table>]");
107 GMT_Usage (API, -2, "Dump level information to <table> or standard output if not given.");
108 gmt_outgrid_syntax (API, 'G', "Create an equalized output grid file called <outgrid>");
109 GMT_Usage (API, 1, "\n-N[<norm>]");
110 GMT_Usage (API, -2, "Use with -G to make an output grid file with standard normal scores. "
111 "Alternatively, append <norm> to normalize the scores to -<norm>/+<norm>.");
112 GMT_Usage (API, 1, "\n-Q Use quadratic equalization scaling [Default is linear].");
113 GMT_Option (API, "R,V,h,.");
114
115 return (GMT_MODULE_USAGE);
116 }
117
parse(struct GMT_CTRL * GMT,struct GRDHISTEQ_CTRL * Ctrl,struct GMT_OPTION * options)118 static int parse (struct GMT_CTRL *GMT, struct GRDHISTEQ_CTRL *Ctrl, struct GMT_OPTION *options) {
119 /* This parses the options provided to grdhisteq and sets parameters in Ctrl.
120 * Note Ctrl has already been initialized and non-zero default values set.
121 * Any GMT common options will override values set previously by other commands.
122 * It also replaces any file names specified as input or output with the data ID
123 * returned when registering these sources/destinations with the API.
124 */
125
126 unsigned int n_errors = 0, n_files = 0;
127 int sval;
128 struct GMT_OPTION *opt = NULL;
129 struct GMTAPI_CTRL *API = GMT->parent;
130
131 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
132
133 switch (opt->option) {
134 case '<': /* Input file (only one is accepted) */
135 if (n_files++ > 0) {n_errors++; continue; }
136 Ctrl->In.active = true;
137 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
138 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
139 break;
140
141 /* Processes program-specific parameters */
142
143 case 'C': /* Get # of cells */
144 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
145 Ctrl->C.active = true;
146 sval = atoi (opt->arg);
147 n_errors += gmt_M_check_condition (GMT, sval <= 0, "Option -C: n_cells must be positive\n");
148 Ctrl->C.value = sval;
149 break;
150 case 'D': /* Dump info to file or stdout */
151 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
152 Ctrl->D.active = true;
153 if (opt->arg[0]) Ctrl->D.file = strdup (opt->arg);
154 break;
155 case 'G': /* Output file for equalized grid */
156 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
157 Ctrl->G.active = true;
158 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
159 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
160 break;
161 case 'N': /* Get normalized scores */
162 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
163 Ctrl->N.active = true;
164 if (opt->arg[0]) Ctrl->N.norm = atof (opt->arg);
165 break;
166 case 'Q': /* Use quadratic scaling */
167 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
168 Ctrl->Q.active = true;
169 break;
170
171 default: /* Report bad options */
172 n_errors += gmt_default_error (GMT, opt->option);
173 break;
174 }
175 }
176
177 n_errors += gmt_M_check_condition (GMT, n_files > 1, "Must specify a single input grid file\n");
178 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input grid file\n");
179 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->G.active, "Option -N: Must also specify output grid file with -G\n");
180 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->Q.active, "Option -N: Cannot be combined with -Q\n");
181 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->C.value <= 0, "Option -C: n_cells must be positive\n");
182 n_errors += gmt_M_check_condition (GMT, !Ctrl->D.active && !Ctrl->G.active, "Either -D or -G is required for output\n");
183 n_errors += gmt_M_check_condition (GMT, Ctrl->In.file && !strcmp (Ctrl->In.file, "="), "Piping of input grid file not supported!\n");
184
185 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
186 }
187
grdhisteq_get_cell(gmt_grdfloat x,struct GRDHISTEQ_CELL * cell,unsigned int n_cells_m1,unsigned int last_cell)188 GMT_LOCAL gmt_grdfloat grdhisteq_get_cell (gmt_grdfloat x, struct GRDHISTEQ_CELL *cell, unsigned int n_cells_m1, unsigned int last_cell) {
189 unsigned int low, high, i;
190 /* Quick bisector search for relevant cell */
191 low = 0;
192 high = n_cells_m1;
193 i = last_cell;
194
195 do {
196 if (cell[i].low <= x && cell[i].high >= x) {
197 return ((gmt_grdfloat)i);
198 }
199 else if (cell[low].low <= x && cell[low].high >= x) {
200 return ((gmt_grdfloat)low);
201 }
202 else if (cell[high].low <= x && cell[high].high >= x) {
203 return ((gmt_grdfloat)high);
204 }
205 else if (cell[i].low > x) {
206 high = i;
207 i = (low + high) / 2;
208 }
209 else if (cell[i].high < x) {
210 low = i;
211 i = (low + high) / 2;
212 }
213 } while (true);
214 return (0.0); /* Cannot get here - just used to quiet compiler */
215 }
216
grdhisteq_do_hist_equalization_cart(struct GMT_CTRL * GMT,struct GMT_GRID * Grid,char * outfile,unsigned int n_cells,bool quadratic)217 GMT_LOCAL struct GRDHISTEQ_CELL *grdhisteq_do_hist_equalization_cart (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, char *outfile, unsigned int n_cells, bool quadratic) {
218 /* Do basic Cartesian histogram equalization */
219 uint64_t i, j, nxy;
220 unsigned int n_cells_m1 = 0, current_cell, pad[4];
221 double delta_cell, target;
222 struct GRDHISTEQ_CELL *cell = NULL;
223 struct GMT_GRID *Orig = NULL;
224
225 cell = gmt_M_memory (GMT, NULL, n_cells, struct GRDHISTEQ_CELL);
226
227 /* Sort the data and find the division points */
228
229 gmt_M_memcpy (pad, Grid->header->pad, 4, int); /* Save the original pad */
230 gmt_grd_pad_off (GMT, Grid); /* Undo pad if one existed so we can sort the entire grid */
231 if (outfile && (Orig = GMT_Duplicate_Data (GMT->parent, GMT_IS_GRID, GMT_DUPLICATE_DATA, Grid)) == NULL) { /* Must keep original if readonly */
232 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Grid duplication failed - memory error?\n");
233 gmt_M_free (GMT, cell);
234 return (NULL);
235 }
236 gmt_sort_array (GMT, Grid->data, Grid->header->nm, GMT_FLOAT);
237
238 nxy = Grid->header->nm;
239 while (nxy > 0 && gmt_M_is_fnan (Grid->data[nxy-1])) nxy--; /* Only deal with real numbers since NaNs will be at end */
240
241 n_cells_m1 = n_cells - 1;
242 current_cell = 0;
243 i = 0;
244 delta_cell = ((double)nxy) / ((double)n_cells);
245
246 while (current_cell < n_cells) {
247
248 if (current_cell == n_cells_m1)
249 j = nxy - 1;
250 else if (quadratic) { /* Use y = 2x - x**2 scaling */
251 target = (current_cell + 1.0) / n_cells;
252 j = lrint (floor (nxy * (1.0 - sqrt (1.0 - target))));
253 }
254 else /* Use simple linear scale */
255 j = lrint (floor ((current_cell + 1) * delta_cell)) - 1;
256
257 cell[current_cell].low = Grid->data[i];
258 cell[current_cell].high = Grid->data[j];
259 cell[current_cell].row = current_cell;
260
261 i = j;
262 current_cell++;
263 }
264
265 if (outfile) { /* Must re-read the grid and evaluate since it got sorted and trodden on... */
266 unsigned int last_cell = n_cells / 2;
267 for (i = 0; i < Grid->header->nm; i++) Grid->data[i] = (gmt_M_is_fnan (Orig->data[i])) ? GMT->session.f_NaN : grdhisteq_get_cell (Orig->data[i], cell, n_cells_m1, last_cell);
268 if (GMT_Destroy_Data (GMT->parent, &Orig) != GMT_NOERROR) {
269 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to free Orig\n");
270 }
271 }
272
273 gmt_grd_pad_on (GMT, Grid, pad); /* Reinstate the original pad */
274
275 return (cell); /* Pass out the cell boundaries */
276 }
277
grdhisteq_do_hist_equalization_geo(struct GMT_CTRL * GMT,struct GMT_GRID * Grid,char * outfile,unsigned int n_cells,bool quadratic)278 GMT_LOCAL struct GRDHISTEQ_CELL *grdhisteq_do_hist_equalization_geo (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, char *outfile, unsigned int n_cells, bool quadratic) {
279 /* Do basic area-weighted histogram equalization */
280 uint64_t i, j, node, nxy = 0;
281 unsigned int n_cells_m1 = 0, current_cell, row, col;
282 double cell_w, delta_w, target_w, wsum = 0.0;
283 struct GRDHISTEQ_CELL *cell = gmt_M_memory (GMT, NULL, n_cells, struct GRDHISTEQ_CELL);
284 struct GMT_GRID *W = gmt_duplicate_grid (GMT, Grid, GMT_DUPLICATE_ALLOC);
285 struct GMT_OBSERVATION *pair = gmt_M_memory (GMT, NULL, Grid->header->nm, struct GMT_OBSERVATION);
286
287 /* Determine the area weights */
288 gmt_get_cellarea (GMT, W);
289 /* Fill in the observation (w,z) observation pairs and find # of points nxy */
290 gmt_M_grd_loop (GMT, Grid, row, col, node) {
291 if (gmt_M_is_fnan (Grid->data[node]) || gmt_M_is_dnan (W->data[node]))
292 continue;
293 pair[nxy].value = Grid->data[node];
294 pair[nxy++].weight = W->data[node];
295 wsum += W->data[node];
296 }
297 gmt_free_grid (GMT, &W, true); /* Done with the area weights grid */
298 /* Sort observations on z */
299 qsort (pair, nxy, sizeof (struct GMT_OBSERVATION), gmtlib_compare_observation);
300 /* Compute normalized cumulative weights */
301 wsum = 1.0 / wsum; /* Do avoid division later */
302 pair[0].weight *= (gmt_grdfloat)wsum;
303 for (i = 1; i < nxy; i++) {
304 pair[i].weight *= (gmt_grdfloat)wsum;
305 pair[i].weight += pair[i-1].weight;
306 }
307
308 /* Find the division points using the normalized 0-1 weights */
309
310 n_cells_m1 = n_cells - 1;
311 current_cell = 0;
312 i = j = 0;
313 cell_w = delta_w = 1.0 / n_cells;
314 while (current_cell < n_cells) {
315
316 if (current_cell == n_cells_m1)
317 j = nxy - 1; /* End at the last sorted point */
318 else if (quadratic) { /* Use y = 2x - x**2 scaling to stretch the target weight at end of box */
319 target_w = (1.0 - sqrt (1.0 - cell_w));
320 while (j < nxy && pair[j].weight < target_w) j++;
321 }
322 else /* Use simple linear scale */
323 while (j < nxy && pair[j].weight < cell_w) j++;
324
325 cell[current_cell].low = pair[i].value;
326 cell[current_cell].high = pair[j].value;
327 cell[current_cell].row = current_cell;
328
329 i = j;
330 current_cell++;
331 cell_w += delta_w;
332 }
333
334 if (outfile) { /* Evaluate grid given its original values */
335 unsigned int last_cell = n_cells / 2;
336 gmt_M_grd_loop (GMT, Grid, row, col, node) {
337 Grid->data[node] = (gmt_M_is_fnan (Grid->data[node])) ? GMT->session.f_NaN : grdhisteq_get_cell (Grid->data[node], cell, n_cells_m1, last_cell);
338 }
339 }
340
341 gmt_M_free (GMT, pair);
342
343 return (cell); /* Pass out the cell boundaries */
344 }
345
grdhisteq_do_hist_equalization(struct GMT_CTRL * GMT,struct GMT_GRID * Grid,char * outfile,unsigned int n_cells,bool quadratic)346 GMT_LOCAL struct GRDHISTEQ_CELL * grdhisteq_do_hist_equalization (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, char *outfile, unsigned int n_cells, bool quadratic) {
347 struct GRDHISTEQ_CELL *C;
348 if (gmt_M_is_geographic (GMT, GMT_IN))
349 C = grdhisteq_do_hist_equalization_geo (GMT, Grid, outfile, n_cells, quadratic);
350 else
351 C = grdhisteq_do_hist_equalization_cart (GMT, Grid, outfile, n_cells, quadratic);
352 return (C);
353 }
354
grdhisteq_compare_indexed_floats(const void * point_1,const void * point_2)355 GMT_LOCAL int grdhisteq_compare_indexed_floats (const void *point_1, const void *point_2) {
356 if (((struct INDEXED_DATA *)point_1)->x < ((struct INDEXED_DATA *)point_2)->x) return (-1);
357 if (((struct INDEXED_DATA *)point_1)->x > ((struct INDEXED_DATA *)point_2)->x) return (+1);
358 return (0);
359 }
360
grdhisteq_compare_indices(const void * point_1,const void * point_2)361 GMT_LOCAL int grdhisteq_compare_indices (const void *point_1, const void *point_2) {
362 if (((struct INDEXED_DATA *)point_1)->i < ((struct INDEXED_DATA *)point_2)->i) return (-1);
363 if (((struct INDEXED_DATA *)point_1)->i > ((struct INDEXED_DATA *)point_2)->i) return (1);
364 return (0);
365 }
366
grdhisteq_do_gaussian_scores(struct GMT_CTRL * GMT,struct GMT_GRID * Grid,double norm)367 GMT_LOCAL int grdhisteq_do_gaussian_scores (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, double norm) {
368 /* Make an output grid file with standard normal scores */
369 unsigned int row, col;
370 uint64_t i = 0, j = 0, ij, nxy;
371 double dnxy;
372 struct INDEXED_DATA *indexed_data = NULL;
373
374 indexed_data = gmt_M_memory (GMT, NULL, Grid->header->nm, struct INDEXED_DATA);
375
376 nxy = Grid->header->nm;
377 gmt_M_grd_loop (GMT, Grid, row, col, ij) {
378 if (gmt_M_is_fnan (Grid->data[ij])) { /* Put NaNs in the back */
379 nxy--;
380 indexed_data[nxy].i = ij;
381 indexed_data[nxy].x = Grid->data[ij];
382 }
383 else {
384 indexed_data[j].i = ij;
385 indexed_data[j].x = Grid->data[ij];
386 j++;
387 }
388 }
389
390 /* Sort on data value */
391
392 qsort (indexed_data, nxy, sizeof (struct INDEXED_DATA), grdhisteq_compare_indexed_floats);
393
394 dnxy = 1.0 / (nxy + 1);
395
396 if (norm != 0.0) norm /= fabs (gmt_zcrit (GMT, dnxy)); /* Normalize by abs(max score) */
397
398 for (i = 0; i < nxy; i++) {
399 indexed_data[i].x = (gmt_grdfloat)gmt_zcrit (GMT, (i + 1.0) * dnxy);
400 if (norm != 0.0) indexed_data[i].x *= (gmt_grdfloat)norm;
401 }
402
403 /* Sort on data index */
404
405 qsort (indexed_data, Grid->header->nm, sizeof (struct INDEXED_DATA), grdhisteq_compare_indices);
406
407 i = 0;
408 gmt_M_grd_loop (GMT, Grid, row, col, ij) Grid->data[ij] = indexed_data[i++].x; /* Load up the grid */
409
410 gmt_M_free (GMT, indexed_data);
411 return (0);
412 }
413
414 #define bailout(code) {gmt_M_free_options (mode); return (code);}
415 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
416
GMT_grdhisteq(void * V_API,int mode,void * args)417 EXTERN_MSC int GMT_grdhisteq (void *V_API, int mode, void *args) {
418 int error = 0;
419
420 double wesn[4];
421
422 struct GMT_GRID *Grid = NULL, *Out = NULL;
423 struct GRDHISTEQ_CELL *Cell = NULL;
424 struct GRDHISTEQ_CTRL *Ctrl = NULL;
425 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
426 struct GMT_OPTION *options = NULL;
427 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
428
429 /*----------------------- Standard module initialization and parsing ----------------------*/
430
431 if (API == NULL) return (GMT_NOT_A_SESSION);
432 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
433 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
434
435 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
436
437 /* Parse the command-line arguments */
438
439 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 */
440 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
441 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
442 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
443
444 /*---------------------------- This is the grdhisteq main code ----------------------------*/
445
446 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid\n");
447 gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Current -R setting, if any */
448 if ((Grid = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) {
449 Return (API->error);
450 }
451 if (gmt_M_is_subset (GMT, Grid->header, wesn)) { /* Subset requested; make sure wesn matches header spacing */
452 if ((error = gmt_M_err_fail (GMT, gmt_adjust_loose_wesn (GMT, wesn, Grid->header), "")))
453 Return (error);
454 }
455 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file, Grid) == NULL) { /* Get subset */
456 Return (API->error);
457 }
458 (void)gmt_set_outgrid (GMT, Ctrl->In.file, false, 0, Grid, &Out); /* true if input is a read-only array */
459 gmt_grd_init (GMT, Out->header, options, true);
460
461 if (Ctrl->N.active)
462 grdhisteq_do_gaussian_scores (GMT, Out, Ctrl->N.norm);
463 else { /* Do histogram equalization and return the cell values via structure array Cell */
464 if ((Cell = grdhisteq_do_hist_equalization (GMT, Out, Ctrl->G.file, Ctrl->C.value, Ctrl->Q.active)) == NULL)
465 Return (GMT_RUNTIME_ERROR); /* Some error */
466 }
467 if (Ctrl->G.active) { /* Output the equalized grid */
468 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) Return (API->error);
469 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Out) != GMT_NOERROR) {
470 Return (API->error);
471 }
472 }
473
474 if (Ctrl->D.active) { /* Wants to dump the equal area cell boundaries */
475 struct GMT_DATASET *D = NULL;
476 struct GMT_DATASEGMENT *S = NULL;
477 uint64_t dim[GMT_DIM_SIZE] = {1, 1, Ctrl->C.value, 3}, row;
478 if ((D = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
479 Return (API->error);
480 }
481 S = D->table[0]->segment[0]; /* Short-hand to the single segment */
482 for (row = 0; row < S->n_rows; row++) {
483 S->data[GMT_X][row] = Cell[row].low;
484 S->data[GMT_Y][row] = Cell[row].high;
485 S->data[GMT_Z][row] = Cell[row].row;;
486 }
487 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_NORMAL, NULL, Ctrl->D.file, D) != GMT_NOERROR) {
488 Return (API->error);
489 }
490 }
491
492 gmt_M_free (GMT, Cell);
493
494 Return (GMT_NOERROR);
495 }
496