grdfill.c

/*--------------------------------------------------------------------
 *
 *	Copyright (c) 1991-2025 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
 *--------------------------------------------------------------------*/
/*
 * Brief synopsis: grdfill.c reads a grid that has one or more holes in
 * it flagged by NaNs and interpolates across the holes given a choice
 * of algorithm.
 *
 * Author:	Paul Wessel
 * Date:	7-JAN-2018
 * Version:	6 API
 */

#include "gmt_dev.h"
#include "longopt/grdfill_inc.h"

#define THIS_MODULE_CLASSIC_NAME	"grdfill"
#define THIS_MODULE_MODERN_NAME	"grdfill"
#define THIS_MODULE_LIB		"core"
#define THIS_MODULE_PURPOSE	"Interpolate across holes in a grid"
#define THIS_MODULE_KEYS	"<G{,>?}"
#define THIS_MODULE_NEEDS	""
#define THIS_MODULE_OPTIONS	"-RVf"

enum GRDFILL_mode {
	ALG_CONSTANT,
	ALG_NN,
	ALG_SPLINE,
	ALG_GRID
};

struct GRDFILL_CTRL {
	struct GRDFILL_In {
		bool active;
		char *file;
	} In;
	struct GRDFILL_A {	/* -Ac|g|n|s[<options>] */
		bool active;
		unsigned int mode;
		double value;
		char *file;
	} A;
	struct GRDFILL_G {	/* -G<outgrid> */
		bool active;
		char *file;
	} G;
	struct GRDFILL_L {	/* -L[p] */
		bool active;
		unsigned int mode;	/* 0 = region, 1 = polygons */
	} L;
	struct GRDFILL_N {	/* -N<value> */
		bool active;
		float value;
	} N;
};

static void *New_Ctrl (struct GMT_CTRL *GMT) {	/* Allocate and initialize a new control structure */
	struct GRDFILL_CTRL *C;

	C = gmt_M_memory (GMT, NULL, 1, struct GRDFILL_CTRL);

	/* Initialize values whose defaults are not 0/false/NULL */

	return (C);
}

static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDFILL_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->G.file);
	gmt_M_free (GMT, C);
}

static int usage (struct GMTAPI_CTRL *API, int level) {
	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
	GMT_Usage (API, 0, "usage: %s %s -Ac|g|n|s[<arg>] | -L[p] [-G%s] [-N<value>] [%s] [%s] [%s] [%s]\n",
		name, GMT_INGRID, GMT_OUTGRID, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_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 grid with NaN holes");
	GMT_Usage (API, 1, "\n-Ac|n|s[<arg>]");
	GMT_Usage (API, -2, "Specify algorithm and any required parameters for in-fill:");
	GMT_Usage (API, 3, "c: Fill in NaN holes with the given constant <value>.");
	GMT_Usage (API, 3, "g: Fill in NaN holes by sampling the gridfile whose name is appended.");
	GMT_Usage (API, 3, "n: Fill in NaN holes with nearest neighbor values; "
		"append <max_radius> nodes for the outward search "
		"[Default radius is sqrt(nx^2+ny^2), with (nx,ny) the dimensions of the grid].");
	GMT_Usage (API, 3, "s: Fill in NaN holes with a spline (optionally append tension).");
	GMT_Usage (API, -2, "Note: -A is required unless -L is used.");
	GMT_Message (API, GMT_TIME_NONE, "\n  OPTIONAL ARGUMENTS:\n");
	gmt_outgrid_syntax (API, 'G', "Give filename for where to write the output grid");
	GMT_Usage (API, 1, "\n-L[p]");
	GMT_Usage (API, -2, "Just list the sub-regions w/e/s/n of each hole. "
		"No grid fill takes place and -G is ignored. "
		"Optionally, append p to write polygons corresponding to these regions.");
	GMT_Usage (API, 1, "\n-N<value>");
	GMT_Usage (API, -2, "Set alternate node <value> to indicate a hole [Default looks for NaN-nodes].");
	GMT_Option (API, "R,V,f,.");

	return (GMT_MODULE_USAGE);
}

static int parse (struct GMT_CTRL *GMT, struct GRDFILL_CTRL *Ctrl, struct GMT_OPTION *options) {
	/* This parses the options provided to grdfill and sets parameters in CTRL.
	 * 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;
	struct GMT_OPTION *opt = NULL;
	struct GMTAPI_CTRL *API = GMT->parent;

	for (opt = options; opt; opt = opt->next) {
		switch (opt->option) {

			case '<':	/* Input and Output files */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->In.active);
				n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file));
				break;
			case '>':	/* Output file  */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
				n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file));
				break;

			/* Processes program-specific parameters */

			case 'A':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
				switch (opt->arg[0]) {
					case 'c':	/* Constant fill */
						Ctrl->A.mode = ALG_CONSTANT;
						Ctrl->A.value = atof (&opt->arg[1]);
						break;
					case 'g':	/* Sample a (coarser?) grid */
						Ctrl->A.mode = ALG_GRID;
						n_errors += gmt_get_required_file (GMT, &opt->arg[1], opt->option, 0, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->A.file));
						break;
					case 'n':	/* Nearest neighbor fill */
						Ctrl->A.mode = ALG_NN;
						Ctrl->A.value = (opt->arg[1]) ? atof (&opt->arg[1]) : -1;
						break;
					case 's':	/* Spline fill */
						Ctrl->A.mode = ALG_SPLINE;
						if (opt->arg[1]) Ctrl->A.value =  atof (&opt->arg[1]);
						break;
					default:
						GMT_Report (API, GMT_MSG_ERROR, "Option -A: Unrecognized algorithm (%c)\n", opt->arg[0]);
						n_errors++;
				}
				break;

			case 'G':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
				if (Ctrl->G.file) {
					GMT_Report (API, GMT_MSG_ERROR, "Specify only one output file\n");
					n_errors++;
				}
				else
					n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file));
				break;

			case 'L':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
				if (opt->arg[0] == 'p')
					Ctrl->L.mode = 1;
				break;

			case 'N':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
				if (opt->arg[0] && !strchr ("Nn", opt->arg[0]))
					Ctrl->N.value = (float) atof (opt->arg);
				else {
					GMT_Report (API, GMT_MSG_ERROR, "Option -N: No value (or NaN) given\n");
					n_errors++;
				}
				break;

			default:	/* Report bad options */
				n_errors += gmt_default_option_error (GMT, opt);
				break;
		}
	}

	n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input grid file\n");
	n_errors += gmt_M_check_condition (GMT, !(Ctrl->A.active || Ctrl->L.active), "Must specify an algorithm with -A unless -L is used\n");
	n_errors += gmt_M_check_condition (GMT, !(Ctrl->L.active || Ctrl->G.file), "Must specify output grid file\n");
	n_errors += gmt_M_check_condition (GMT, !(Ctrl->A.active || Ctrl->L.active), "Must specify either -A or -L\n");

	return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}

GMT_LOCAL int grdfill_do_constant_fill (struct GMT_GRID *G, unsigned int limit[], gmt_grdfloat value) {
	/* Algorithm 1: Replace NaNs with a constant value */
	uint64_t node;

	for (unsigned int row = limit[YLO]; row <= limit[YHI]; row++) {
		for (unsigned int col = limit[XLO]; col <= limit[XHI]; col++) {
			node = gmt_M_ijp (G->header, row, col);
			if (gmt_M_is_fnan (G->data[node]))
				G->data[node] = value;
		}
	}
	return GMT_NOERROR;
}

#if 1
GMT_LOCAL int grdfill_do_splinefill (struct GMTAPI_CTRL *API, struct GMT_GRID *G, double wesn[], unsigned int limit[], unsigned int n_in_hole, double value) {
	/* Algorithm 2: Replace NaNs with a spline */
	char input[GMT_VF_LEN] = {""}, output[GMT_VF_LEN] = {""}, args[GMT_LEN256] = {""}, method[GMT_LEN32] = {""};
	unsigned int row, col, row_hole, col_hole, mode, d_limit[4], n_constraints;
	uint64_t node, node_hole, k = 0, dim[GMT_DIM_SIZE] = {0, 0, 0, 0};
	double *x = NULL, *y = NULL;
	gmt_grdfloat *z = NULL;
	struct GMT_VECTOR *V = NULL;
	struct GMT_GRID *G_hole = NULL;
	struct GMT_CTRL *GMT = API->GMT;

	/* Allocate a vector container for input to greenspline */
	dim[0] = 3;	/* Want three input columns but let length be 0 - this signals that no vector allocations should take place */
	if ((V = GMT_Create_Data (API, GMT_IS_VECTOR, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
		return (API->error);
	}
	/* Create a virtual file to hold the resampled grid */
	if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output) == GMT_NOTSET) {
		return (API->error);
	}
	/* Add up to 2 rows/cols around hole, but watch for grid edges */
	gmt_M_memcpy (d_limit, limit, 4, unsigned int);	/* d_limit will be used to set the grid domain */
	/* Undo the half grid inc padding done in the main */
	wesn[XLO] += 0.5 * G->header->inc[GMT_X];
	wesn[XHI] -= 0.5 * G->header->inc[GMT_X];
	wesn[YLO] += 0.5 * G->header->inc[GMT_Y];
	wesn[YHI] -= 0.5 * G->header->inc[GMT_Y];
	for (k = 1; k <= 2; k++) {
		if (d_limit[XLO]) d_limit[XLO]--, wesn[XLO] -= G->header->inc[GMT_X];	/* Move one column westward */
		if (d_limit[XHI] < (G->header->n_columns-1)) d_limit[XHI]++, wesn[XHI] += G->header->inc[GMT_X];	/* Move one column eastward */
		if (d_limit[YLO]) d_limit[YLO]--, wesn[YHI] += G->header->inc[GMT_Y];	/* Move one row northward */
		if (d_limit[YHI] < (G->header->n_rows-1)) d_limit[YHI]++, wesn[YLO] -= G->header->inc[GMT_Y];	/* Move one row southward */
	}
	n_constraints = (d_limit[YHI] - d_limit[YLO] + 1) * (d_limit[XHI] - d_limit[XLO] + 1) - n_in_hole;
	x = gmt_M_memory (GMT, NULL, n_constraints, double);
	y = gmt_M_memory (GMT, NULL, n_constraints, double);
	z = gmt_M_memory (GMT, NULL, n_constraints, gmt_grdfloat);
	for (row = d_limit[YLO], k = 0; row <= d_limit[YHI]; row++) {
		for (col = d_limit[XLO]; col <= d_limit[XHI]; col++) {
			node = gmt_M_ijp (G->header, row, col);
			if (gmt_M_is_fnan (G->data[node])) continue;
			x[k] = gmt_M_grd_col_to_x (GMT, col, G->header);
			y[k] = gmt_M_grd_row_to_y (GMT, row, G->header);
			z[k] = G->data[node];
			k++;
		}
	}
	V->n_rows = n_constraints;	/* Must specify how many input points we have */
	GMT_Put_Vector (API, V, GMT_X, GMT_DOUBLE, x);
	GMT_Put_Vector (API, V, GMT_Y, GMT_DOUBLE, y);
	GMT_Put_Vector (API, V, GMT_Z, GMT_FLOAT,  z);
	/* Associate our input data vectors with a virtual input file */
	if (GMT_Open_VirtualFile (API, GMT_IS_DATASET|GMT_VIA_VECTOR, GMT_IS_POINT, GMT_IN|GMT_IS_REFERENCE, V, input) == GMT_NOTSET)
		return (API->error);
	/* Prepare the greenspline command-line arguments */
	mode = (gmt_M_is_geographic (GMT, GMT_IN)) ? 2 : 1;
	if (value > 0.0)
		sprintf (method, "t%g", value);
	else
		sprintf (method, "c");
	sprintf (args, "%s -G%s -S%s -R%.16g/%.16g/%.16g/%.16g -I%.16g/%.16g -D%d", input, output, method, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], G->header->inc[GMT_X], G->header->inc[GMT_Y], mode);
	if (G->header->registration == GMT_GRID_PIXEL_REG) strcat (args, " -r");
	strcat (args, " --GMT_HISTORY=readonly");
   	/* Run the greenspline module */
	GMT_Report (API, GMT_MSG_INFORMATION, "Calling greenspline with args %s\n", args);
  	if (GMT_Call_Module (API, "greenspline", GMT_MODULE_CMD, args)) {
		return (API->error);	/* Some sort of failure */
	}
	if ((G_hole = GMT_Read_VirtualFile (API, output)) == NULL) {	/* Load in the resampled grid */
		return (API->error);	/* Some sort of failure */
	}
	/* Use 0-nx,0-ny for the hole index and the large grid G nodes for the holes */
	for (row_hole = 0, row = d_limit[YLO]; row_hole < G_hole->header->n_rows; row_hole++, row++) {
		for (col_hole = 0, col = d_limit[XLO]; col_hole < G_hole->header->n_columns; col_hole++, col++) {
			node = gmt_M_ijp (G->header, row, col);
			if (!gmt_M_is_fnan (G->data[node])) continue;	/* Had data value */
			node_hole = gmt_M_ijp (G_hole->header, row_hole, col_hole);
			G->data[node] = G_hole->data[node_hole];	/* Replace the NaN with splined value */
		}
	}

	/* Close the two virtual files */
	if (GMT_Close_VirtualFile (API, output)) {
		GMT_Report (API, GMT_MSG_ERROR, "Failed to close virtual output file %s\n", output);
		return (API->error);
	}
	if (GMT_Close_VirtualFile (API, input)) {
		GMT_Report (API, GMT_MSG_ERROR, "Failed to close virtual input file %s\n", input);
		return (API->error);
	}

	/* Free our custom vectors */
	gmt_M_free (GMT, x);	gmt_M_free (GMT, y);	gmt_M_free (GMT, z);
	if (GMT_Destroy_Data (API, &G_hole) != GMT_NOERROR) {
		GMT_Report (API, GMT_MSG_ERROR, "Failed to destroy temporary hold grid\n");
		return (API->error);
	}
	return (GMT_NOERROR);
}
#endif

GMT_LOCAL unsigned int grdfill_trace_the_hole (struct GMT_GRID *G, uint64_t node, unsigned int row, unsigned int col, int64_t *step, char *ID, unsigned int limit[]) {
	/* Determine all the direct neighbor nodes in the W/E/S/N directions that are also NaN, recursively.
	 * This is a limited form of Moore neighborhood called Von Neumann neighborhood since we only do the
	 * four cardinal directions and ignore the diagonals.  Ignoring the diagonal means two holes that
	 * touch at a diagonal will be encoded as two separate holes whereas a full Moore neighborhood
	 * calculation would find a single hole.  For our purposes (filling in the hole), smaller holes
	 * are less computationally intensive, hence we limit ourselves to good old Von Neumann. */
	static int drow[4] = {1, 0, -1, 0}, dcol[4] = {0, 1, 0, -1};	/* Change in row,col per cardinal direction */
	unsigned int side, next_row, next_rcol, n_nodes = 0;
	int64_t ij;

	for (side = 0; side < 4; side++) {	/* For each of the 4 cardinal directions */
		ij = node + step[side];		/* This is the grid node (and ID node) of this cardinal point */
		if (ID[ij] == 0 && gmt_M_is_fnan (G->data[ij])) {	/* Hole extends in this direction, follow it to its next neighbor */
			next_row  = (unsigned int)(row + drow[side]);	/* Set col and row of next point */
			next_rcol = (unsigned int)(col + dcol[side]);
			ID[ij] = 1;	/* Mark this node as part of the current hole */
			if (next_rcol < limit[XLO]) limit[XLO] = next_rcol;	/* Update min/max row and col values */
			else if (next_rcol > limit[XHI]) limit[XHI] = next_rcol;
			if (next_row < limit[YLO]) limit[YLO] = next_row;
			else if (next_row > limit[YHI]) limit[YHI] = next_row;
			/* Recursively trace this nodes next neighbors as well */
			n_nodes = n_nodes + 1 + grdfill_trace_the_hole (G, ij, next_row, next_rcol, step, ID, limit);
		}
	}
	return (n_nodes);
}

GMT_LOCAL int64_t grdfill_find_nearest (int64_t i, int64_t j, int64_t *r2, int64_t *is, int64_t *js, int64_t *xs, int64_t *ys) {
	/* function to find the nearest point based on previous search, smallest distance outside a radius */
	int64_t ct = 0, nx, ny, nx1, ii, k = 0, rr, nx1_2, ny_2;

	rr = INTMAX_MAX;	/* Ensure we reset this the first time */

	/* starting with nx = ny */
	nx = (int64_t)(sqrt((double)(*r2) / 2.0));
	/* loop over possible nx, find smallest rr */
	for (nx1 = nx; nx1 <= (int64_t)sqrt((double)(*r2)) + 1; nx1++) {
		if ((nx1_2 = (nx1 * nx1)) < (*r2))
			ny = (int64_t)(sqrt((double)((*r2) - nx1_2)));
		else
			ny = 0;
		while ((nx1_2 + ny * ny ) <= (*r2) && ny <= nx1)
			ny++;
		ny_2 = ny * ny;
		if (ny <= nx1) {
			if (rr > (nx1_2 + ny_2)) {
				k = 0;
				rr = nx1_2 + ny_2;
				xs[0] = nx1;
				ys[0] = ny;
			}
			else if (rr == (nx1_2 + ny_2)) {
				k++;
				xs[k] = nx1;
				ys[k] = ny;
			}
		}
	}

	/* Return the grid index, changing the order may lead to different solution, but mainly because nearest-neighbor is non-unique */
	for (ii = 0; ii <= k; ii++) {
		nx = xs[ii];
		ny = ys[ii];

		if (ny == 0) {
			js[ct] =   0;	is[ct++] =  nx;
			js[ct] =   0;	is[ct++] = -nx;
			js[ct] =  nx;	is[ct++] =   0;
			js[ct] = -nx;	is[ct++] =   0;
		}
		else if (nx != ny){
			js[ct] =  ny;	is[ct++] =  nx;
			js[ct] =  ny;	is[ct++] = -nx;
			js[ct] = -ny;	is[ct++] =  nx;
			js[ct] = -ny;	is[ct++] = -nx;
			js[ct] =  nx;	is[ct++] =  ny;
			js[ct] = -nx;	is[ct++] =  ny;
			js[ct] =  nx;	is[ct++] = -ny;
			js[ct] = -nx;	is[ct++] = -ny;
		}
		else {
			js[ct] =  nx;	is[ct++] =  nx;
			js[ct] =  nx;	is[ct++] = -nx;
			js[ct] = -nx;	is[ct++] =  nx;
			js[ct] = -nx;	is[ct++] = -nx;
		}
 	}

	for (ii = 0; ii < ct; ii++) {	/* Set absolute row/col values */
		is[ii] = is[ii] + i;
		js[ii] = js[ii] + j;
	}

	*r2 = rr;	/* Update radius */

	/* return the number of indexes of the nearest neighbor */
	return (ct);
}

GMT_LOCAL void grdfill_nearest_interp (struct GMT_CTRL *GMT, struct GMT_GRID *In, struct GMT_GRID *Out, int64_t radius) {
	uint64_t ij, node;
	int64_t nx = In->header->n_columns, ny = In->header->n_rows;
 	int64_t i, j, flag, ct, k, recx = 1, recy = 1, cs = 0, rr;
 	int64_t *is = NULL, *js = NULL, *xs = NULL, *ys = NULL;
	gmt_grdfloat *m = In->data, *m_interp = Out->data;	/* Short-hand for input and output grids */
	double rad2;

 	/* Allocate memory for temporary indexes */
 	is = gmt_M_memory (GMT, NULL, 2048, int64_t);
 	js = gmt_M_memory (GMT, NULL, 2048, int64_t);
 	xs = gmt_M_memory (GMT, NULL, 512, int64_t);
 	ys = gmt_M_memory (GMT, NULL, 512, int64_t);

	if (radius == -1)	/* Set default radius */
		radius = (int64_t)floor (sqrt ((double)(nx*nx + ny*ny)));
	rad2 = (double)(radius * radius);

	GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Interpolating to nearest neighbor...\n");
	gmt_M_row_loop (GMT, In, i) {	/* Loop over each row in grid */
 		GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Working on row %" PRIi64 "\n", i);
		rr = 0;
		gmt_M_col_loop (GMT, In, i, j, ij) {	/* Loop over all columns */
			if (!gmt_M_is_fnan (m[ij]))	/* Already duplicated in the calling program */
				rr = 0;
			else {	/* search nearest neighbor, use previous nearest distance to exclude certain search area. */
				flag = 0;
				/* set the starting search radius based on last nearest distance */
				if (rr >= 4) {
					if (recy > 0 && recx > 0)
						rr = (recx-1)*(recx-1)+(recy-1)*(recy-1)-1;
					else if (recy == 0 && recx > 0)
						rr = (recx-1)*(recx-1)-1;
					else if (recy > 0 && recx == 0)
						rr = (recy-1)*(recy-1)-1;
				}
				else
					rr = 0;

				while (flag == 0 && rr <= rad2) {
					ct = grdfill_find_nearest (i, j, &rr, is, js, xs, ys);
					cs++;
					for (k = 0; k < ct; k++) {
						if (is[k] >= 0 && is[k] < ny && js[k] >=0 && js[k] < nx) {
							node = gmt_M_ijp (In->header, is[k], js[k]);
							if (!gmt_M_is_fnan (m[node])) {
								m_interp[ij] = m[node];
								flag = 1;
								recx = int64_abs (is[k]-i);
								recy = int64_abs (js[k]-j);
								break;
							}
						}
					}
				}
			}
			if (i == 0 && rr == -1)
				GMT_Report (GMT->parent, GMT_MSG_DEBUG, "(%d %d %d %d)\n", j, rr, recx, recy);
		}
	}

	GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%" PRIi64 " number of searches used\n", cs);

	gmt_M_free (GMT, is);
	gmt_M_free (GMT, js);
	gmt_M_free (GMT, xs);
	gmt_M_free (GMT, ys);
}

GMT_LOCAL int grdfill_sample (struct GMT_CTRL *GMT, struct GMT_GRID *In, struct GMT_GRID *Out, char *file) {
	/* Algorithm 3: Replace NaNs with values sampled from another grid */
	/* 1. Create x,y arrays of NaN locations in the grid and an empty z array
	 * 2. Call grdtrack with virtual input and virtual output
	 * 3. Fill in the holes in the grid with these values
	 */
	openmp_int row, col, k;
	uint64_t dim[4] = {0, 0, 0, 0}, node, n_NaNs = 0, n_bad = 0, n_alloc = 0;
	char input[GMT_VF_LEN] = {""}, output[GMT_VF_LEN] = {""}, args[GMT_LEN256] = {""};
	double *data[3] = {NULL, NULL, NULL};
	struct GMT_VECTOR *V_in = NULL, *V_out = NULL;
	struct GMTAPI_CTRL *API = GMT->parent;	/* Shorthand */

	/* Build array list of node coordinates where we have NaNs */
	gmt_M_grd_loop (GMT, In, row, col, node) {
		if (!gmt_M_is_fnan (In->data[node])) continue;	/* Not part of a hole */
		if (n_NaNs >= n_alloc) {	/* Allocate more memory (including the first time) */
			if (n_alloc == 0) n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; else n_alloc *= 2;
			for (k = 0; k < 2; k++) data[k] = gmt_M_memory (GMT, data[k], n_alloc, double);
		}
		/* Keep the coordinates of the NaN node */
		data[GMT_X][n_NaNs] = In->x[col];
		data[GMT_Y][n_NaNs] = In->y[row];
		n_NaNs++;	/* Number of NaNs so far */
	}
	if (n_NaNs == 0) {	/* Well, that was a waste of our time */
		GMT_Report (API, GMT_MSG_WARNING, "There are no NaNs in your grid - returning the input grid as output\n");
		for (k = 0; k < 2; k++) gmt_M_free (GMT, data[k]);
		return (GMT_NOERROR);
	}
	/* Finalize the memory allocation (and z for the first time) */
	for (k = 0; k < 3; k++) data[k] = gmt_M_memory (GMT, data[k], n_NaNs, double);
	/* Allocate two vector containers for input and output separately */
	dim[0] = 2;	/* Want two input columns but let length be 0 - this signals that no vector allocations should take place */
	if ((V_in = GMT_Create_Data (API, GMT_IS_VECTOR, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) return (API->error);
	dim[0] = 3;	/* Want three output columns [we will share the first two with input] */
	if ((V_out = GMT_Create_Data (API, GMT_IS_VECTOR, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) return (API->error);
	V_in->n_rows = V_out->n_rows = n_NaNs;	/* Must specify how many input points we have */
	/* Hook these up to our containers, reusing x,y as both in and out x/y vectors */
	GMT_Put_Vector (API, V_in,  GMT_X, GMT_DOUBLE, data[GMT_X]);
	GMT_Put_Vector (API, V_in,  GMT_Y, GMT_DOUBLE, data[GMT_Y]);
	GMT_Put_Vector (API, V_out, GMT_X, GMT_DOUBLE, data[GMT_X]);
	GMT_Put_Vector (API, V_out, GMT_Y, GMT_DOUBLE, data[GMT_Y]);
	GMT_Put_Vector (API, V_out, GMT_Z, GMT_DOUBLE, data[GMT_Z]);
	if (GMT_Open_VirtualFile (API, GMT_IS_DATASET|GMT_VIA_VECTOR, GMT_IS_POINT, GMT_IN, V_in, input)) {
		GMT_Report (API, GMT_MSG_ERROR, "Failed to use vectors as a virtual dataset!\n");
		return (API->error);
	}
	if (GMT_Open_VirtualFile (API, GMT_IS_DATASET|GMT_VIA_VECTOR, GMT_IS_POINT, GMT_OUT, V_out, output)) {
		GMT_Report (API, GMT_MSG_ERROR, "Failed to use vectors as a virtual dataset!\n");
		return (API->error);
	}
	sprintf (args, "-G%s -R%.16g/%.16g/%.16g/%.16g %s > %s", file, In->header->wesn[XLO], In->header->wesn[XHI],
		In->header->wesn[YLO], In->header->wesn[YHI], input, output);	/* Build grdtrack command and limit with -R */
	if (GMT_Call_Module (API, "grdtrack", GMT_MODULE_CMD, args)) {
		GMT_Report (API, GMT_MSG_ERROR, "Run-time error from grdtrack\n");
		return (API->error);
	}
	n_NaNs = 0;	/* Rewind and march through the grid again */
	gmt_M_grd_loop (GMT, In, row, col, node) {
		if (!gmt_M_is_fnan (In->data[node])) continue;	/* Not part of a hole */
		Out->data[node] = data[GMT_Z][n_NaNs++];	/* Replace the values at NaN-nodes */
		if (gmt_M_is_fnan (Out->data[node])) n_bad++;	/* Oh-oh, not getting a clean output grid */
	}
	if (n_bad)
		GMT_Report (API, GMT_MSG_WARNING, "There are still %" PRIu64 " NaNs remaining after -Ag sampling!\n", n_bad);

	for (k = 0; k < 3; k++) gmt_M_free (GMT, data[k]);	/* Free the temp memory */
	return (GMT_NOERROR);
}

#define bailout(code) {gmt_M_free_options (mode); return (code);}
#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}

EXTERN_MSC int GMT_grdfill (void *V_API, int mode, void *args) {
	char *ID = NULL, *RG_orig_hist = NULL;
	int error = 0, RG_id = 0;
	openmp_int row, col;
	unsigned int hole_number = 0, limit[4], n_nodes;
	uint64_t node, offset;
	int64_t off[4];
	double wesn[4];
	struct GMT_GRID *Grid = NULL;
	struct GMT_RECORD *Out = NULL;
	struct GRDFILL_CTRL *Ctrl = NULL;
	struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
	struct GMT_OPTION *options = NULL;
	struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API);	/* Cast from void to GMTAPI_CTRL pointer */

	/*----------------------- Standard module initialization and parsing ----------------------*/

	if (API == NULL) return (GMT_NOT_A_SESSION);
	if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE));	/* Return the purpose of program */
	options = GMT_Create_Options (API, mode, args);	if (API->error) return (API->error);	/* Set or get option list */

	if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error);	/* Give usage if requested */

	/* Parse the command-line arguments */

	if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, module_kw, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
	if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
	Ctrl = New_Ctrl (GMT);	/* Allocate and initialize a new control structure */
	if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);

	/*---------------------------- This is the grdfill main code ----------------------------*/

	if ((Grid = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY|GMT_GRID_NEEDS_PAD1, NULL, Ctrl->In.file, NULL)) == NULL) {	/* Get header only */
		Return (API->error);
	}

	if (GMT->common.R.active[RSET]) {	/* Specified a subset */
		bool global = false;
		global = gmt_grd_is_global (GMT, Grid->header);
		if (!global && (GMT->common.R.wesn[XLO] < Grid->header->wesn[XLO] || GMT->common.R.wesn[XHI] > Grid->header->wesn[XHI])) error++;
		if (GMT->common.R.wesn[YLO] < Grid->header->wesn[YLO] || GMT->common.R.wesn[YHI] > Grid->header->wesn[YHI]) error++;
		if (error) {
			GMT_Report (API, GMT_MSG_ERROR, "Subset exceeds data domain!\n");
			Return (GMT_RUNTIME_ERROR);
		}
		if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, GMT->common.R.wesn, Ctrl->In.file, Grid) == NULL) {
			Return (API->error);	/* Get subset */
		}
	}
	else if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, Grid) == NULL) {
		Return (API->error);	/* Get all */
	}

	if (Ctrl->N.active) {	/* User wants a specific value to indicate a hole instead of NaN; replace those values with NaNs since that is what is expected below */
		gmt_M_grd_loop (GMT, Grid, row, col, node) {	/* Loop over all grid nodes */
			if (floatAlmostEqualZero (Grid->data[node], Ctrl->N.value))
				Grid->data[node] = GMT->session.f_NaN;
		}
	}
	/* Here any hole is identified as a patch of NaNs */

	if (Ctrl->A.mode == ALG_NN) {	/* Do Eric Xu's NN algorithm and bail */
		int64_t radius = lrint (Ctrl->A.value);
		struct GMT_GRID *New = NULL;
		if ((New = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_DATA, Grid)) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Unable to duplicate input grid!\n");
			Return (API->error);
		}
		grdfill_nearest_interp (GMT, Grid, New, radius);	/* Perform the NN replacements */

		if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, New)) {
			GMT_Report (API, GMT_MSG_ERROR, "Failed to write output grid!\n");
			Return (API->error);
		}
		Return (GMT_NOERROR);
	}

	if (Ctrl->A.mode == ALG_GRID) {	/* Basically a grdtrack virtual file exercise */
		struct GMT_GRID *New = NULL;

		if (gmt_set_outgrid (GMT, Ctrl->In.file, false, 0, Grid, &New))	/* true if input is a read-only array */
			gmt_M_memcpy (New->data, Grid->data, New->header->size, gmt_grdfloat);	/* First duplicate all nodes to the new grid */

		if ((error = grdfill_sample (GMT, Grid, New, Ctrl->A.file)))
			Return (error);	/* Not a happy ending */

		/* Output the filled-in grid */
		if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, New)) {
			GMT_Report (API, GMT_MSG_ERROR, "Failed to write output grid!\n");
			Return (API->error);
		}
		Return (GMT_NOERROR);
	}

	if (Ctrl->L.active) {
		if (GMT_Init_IO (API, GMT_IS_DATASET, (Ctrl->L.mode) ? GMT_IS_POLYGON : GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {	/* Registers default output destination, unless already set */
			Return (API->error);
		}
		if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_OFF) != GMT_NOERROR) {	/* Enables data output and sets access mode */
			Return (API->error);
		}
		if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) {	/* Sets output geometry */
			Return (API->error);
		}
		if ((error = GMT_Set_Columns (API, GMT_OUT, 2 + 2*(1-Ctrl->L.mode), GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
			Return (API->error);
		}
		if (Ctrl->L.mode) gmt_set_segmentheader (GMT, GMT_OUT, true);
	}

	if (Ctrl->A.mode == ALG_SPLINE) {	/* Must preserve the original -RG history which greenspline messes with */
		RG_id = gmt_get_option_id (0, "R") + 1;
		if (GMT->init.history[RG_id]) RG_orig_hist = strdup (GMT->init.history[RG_id]);
	}

	/* To avoid having to check every row,col for being inside the grid we set
	 * the boundary row/cols in the ID grid to 1. */

	ID = gmt_M_memory_aligned (GMT, NULL, Grid->header->size, char);
	/* Set the top and bottom boundary rows to 1 */
	offset = (uint64_t)(Grid->header->pad[YHI] + Grid->header->n_rows) * Grid->header->mx;
	for (node = 0; node < (uint64_t)Grid->header->pad[YHI]*Grid->header->mx; node++) ID[node] = ID[node+offset] = 1;
	/* Set the left and right boundary columns to 1 */
	offset = Grid->header->pad[XLO] + Grid->header->n_columns;
	for (row = 0; row < (openmp_int)Grid->header->my; row++) {
		for (col = 0; col < (openmp_int)Grid->header->pad[XLO]; col++)
			ID[row*Grid->header->mx+col] = 1;
		for (col = 0; col < (openmp_int)Grid->header->pad[XHI]; col++)
			ID[row*Grid->header->mx+offset+col] = 1;
	}
	/* Initiate the node offsets in the cardinal directions */
	off[0] = Grid->header->mx;	off[1] = 1; 	off[2] = -off[0]; off[3] = -off[1];

	Out = gmt_new_record (GMT, wesn, NULL);	/* Since we only need to worry about numerics in this module */
	gmt_M_grd_loop (GMT, Grid, row, col, node) {	/* Loop over all grid nodes */
		if (ID[node]) continue;	/* Already identified as part of a hole */
		if (gmt_M_is_fnan (Grid->data[node])) {	/* Node is part of a new hole */
			limit[XLO] = limit[XHI] = col;	/* Initiate min/max col to this single node */
			limit[YLO] = limit[YHI] = row;	/* Initiate min/max row to this single node */
			ID[node] = 1;	/* Flag this node as part of a hole */
			++hole_number;	/* Increase the current hole number */
			/* Trace all the contiguous neighbors, updating the bounding box as we go along */
			n_nodes = 1 + grdfill_trace_the_hole (Grid, node, row, col, off, ID, limit);
			wesn[XLO] = gmt_M_grd_col_to_x (GMT, limit[XLO], Grid->header) - 0.5 * Grid->header->inc[GMT_X];
			wesn[XHI] = gmt_M_grd_col_to_x (GMT, limit[XHI], Grid->header) + 0.5 * Grid->header->inc[GMT_X];
			wesn[YLO] = gmt_M_grd_row_to_y (GMT, limit[YHI], Grid->header) - 0.5 * Grid->header->inc[GMT_Y];
			wesn[YHI] = gmt_M_grd_row_to_y (GMT, limit[YLO], Grid->header) + 0.5 * Grid->header->inc[GMT_Y];
			GMT_Report (API, GMT_MSG_INFORMATION, "Hole BB %u: -R: %g/%g/%g/%g [%u nodes]\n", hole_number, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], n_nodes);
			if (Ctrl->L.active) {
				if (Ctrl->L.mode) {	/* Write a closed polygon */
					double tmp[4];
					gmt_M_memcpy (tmp, wesn, 4, double);
					GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
					wesn[GMT_X] = tmp[XLO];	wesn[GMT_Y] = tmp[YLO];
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
					wesn[GMT_X] = tmp[XHI];
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
					wesn[GMT_Y] = tmp[YHI];
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
					wesn[GMT_X] = tmp[XLO];
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
					wesn[GMT_Y] = tmp[YLO];
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
				}
				else	/* Write west east south north limits */
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);
			}
			else {
				switch (Ctrl->A.mode) {
					case ALG_CONSTANT:	/* Fill in using a constant value */
						error = grdfill_do_constant_fill (Grid, limit, (gmt_grdfloat)Ctrl->A.value);
						break;
					case ALG_SPLINE:	/* Fill in using a spline */
						error = grdfill_do_splinefill (API, Grid, wesn, limit, n_nodes, Ctrl->A.value);
						break;
				}
				if (error) {
					GMT_Report (API, GMT_MSG_INFORMATION, "Failed to fill hole %u\n", hole_number);
					continue;
				}
			}
		}
	}
	if (hole_number) GMT_Report (API, GMT_MSG_INFORMATION, "Found %u holes\n", hole_number);
	gmt_M_free_aligned (GMT, ID);
	gmt_M_free (GMT, Out);

	if (Ctrl->L.active) {
		if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) {	/* Disables further data output */
			free(RG_orig_hist);
			Return (API->error);
		}
	}
	else {	/* Must write the revised grid whether there are any holes or not */
		if (hole_number == 0) GMT_Report (API, GMT_MSG_WARNING, "No holes detected in grid - grid unchanged\n");
		if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) {
			free(RG_orig_hist);
			Return (API->error);
		}
		if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
			free(RG_orig_hist);
			Return (API->error);
		}
	}

	if (Ctrl->A.mode == ALG_SPLINE) {
		if (GMT->init.history[RG_id]) gmt_M_str_free (GMT->init.history[RG_id]);
		if (RG_orig_hist) GMT->init.history[RG_id] = RG_orig_hist;
	}

	Return (GMT_NOERROR);
}