/* imgbm_polish.c
 *
 * 2007 April 23, WHFS.  Replaces imgbm_polish and includes 
 * constraint ocean/not-ocean bit mask.  Uses old-style
 * img bm constraint pairs.
 * 
 *
 * Reads an input img file, the result of an unconstrained
 * and unpolished gravity prediction of bathymetry.  Creates
 * a new img file which has been "polished" to fit the given
 * imgbm constraints and the given ocean/not-ocean bitmask.
 * The output file is even at unconstrained points and odd at
 * constrained points.
 *
 * The coordinate system must be set to indicate whether we
 * are at 1 or 2 arc-minute and 72 or 81 latitude.
 *
 *
 * The tension parameter is set as usual.  [-T0.6].
 *
 * The periodic boundary conditions are satisfied in the X direction
 * and natural BCs (no bending moment, no shear stress) are used at
 * the N and S edges of the array.  Extra rows are used to hold the
 * N and S BCs.
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
/* #include <float.h>
#include <limits.h>  */

#define MAX_ITER 100	/* Maximum number of iterations to do at any mesh size  */
#define CONVRG 10	/* Meters, must be even integer bigger than 2  */
#define TENSION 0.6	/* tension factor between 0 (biharmonic) and 1 (laplacian)  */
#define SOR 1.4		/* Successive over-relaxation factor  */
#define TALLEST	9000	/* Used to quality-control the constraints  */
#define DEEPEST	-12000	/* Used to quality-control the constraints  */
#define MAX_DELTA 2000	/* Do not use a point at sea if it changes solution by more than this */
#define MOST_SHOAL -6	/* Value not to exceed when iteration would pop ocean point above sea level */

struct	DIMS {
	size_t	nx;
	size_t	ny;
	size_t	nxny;	/* nx * ny  */
	size_t	nxmy;	/* nx * (ny+4); holds 2 top and 2 bottom rows of boundary conditions */
	size_t	nmask;	/* nx * ny / 8; size of bit masks  */
	size_t	ijoff;	/* nx * 2; index to start of data after two rows of BCs  */
} dims;

struct IMGPAIR {
	int     index;
	int     value;
};

static unsigned char maskset[8] = {128, 64, 32, 16, 8, 4, 2, 1};

short int *ain, *a;

unsigned char *hitmask, *seamask;

int	debug = 0;


int	main (int argc, char **argv) {
	
	char	*bfile = NULL, *efile = NULL;
	double	tension = TENSION, f20, f10, f11;
	int	i, error = 0, grid, read_input_img(), read_constraints();
	int	set_dimensions (struct DIMS *dims, int *gcf, int latmax, int minutes);
	int	divide_grid (int grid);
	int	read_landsea_mask (char *fname);
	void	iterate(), initial_reset(), interpolate_new_mesh(), do_output();
	int	start_grid = 0;
	int	latmax = 0;
	int	minutes = 0;
	int	max_iter = MAX_ITER;
	int	converge = CONVRG;
	int	n_tries, last_grid;

	for (i = 1; !error && i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'B':
					bfile = &argv[i][2];
					break;
				case 'C':
					converge = atoi(&argv[i][2]);
					if (converge <= 0 || converge%2) {
						fprintf(stderr,"imgbm_polish:  -C converge must be a positive even integer.\n");
						error++;
					}
					break;
				case 'E':
					efile = &argv[i][2];
					break;
				case 'G':
					start_grid = atoi(&argv[i][2]);
					if (start_grid <= 0) {
						fprintf(stderr,"imgbm_polish:  -G start grid must be positive and a factor of both nx and ny.\n");
						error++;
					}
					break;
				case 'I':
					max_iter = atoi(&argv[i][2]);
					if (max_iter <= 0) {
						fprintf(stderr,"imgbm_polish:  -I max iterations must be a positive integer.\n");
						error++;
					}
					break;
				case 'L':
					latmax = atoi(&argv[i][2]);
					if (latmax != 72 && latmax != 81) {
						fprintf(stderr,"imgbm_polish:  -L latitude maximum must be set to integer 72 or 81.\n");
						error++;
					}
					break;
				case 'M':
					minutes = atoi(&argv[i][2]);
					if (minutes != 1 && minutes != 2) {
						fprintf(stderr,"imgbm_polish:  -M minutes must be set to 1 or 2.\n");
						error++;
					}
					break;
				case 'T':
					tension = atof(&argv[i][2]);
					if (tension < 0.0 || tension > 1.0) {
						fprintf(stderr,"imgbm_polish:  -T tension must be between 0.0 and 1.0.\n");
						error++;
					}
					break;
				default:
					error++;
			}
		}
		else {
			error++;
		}
	}

	if ( (!(efile)) || (!(bfile)) ) error++;
	if (latmax != 72 && latmax != 81) error++;
	if (minutes != 2 && minutes != 1) error++;

	if (error) {
		fprintf(stderr,"usage:  imgbm_polish -B<landsea.bits> -E<input.img> -L<latmax> -M<minutes> [-C<c>] [-G<g>] [-I<i>] [-T<t>] < imgbm_combined_data > output.img\n");
		fprintf (stderr, "-B sets file name for ocean/not-ocean bitmask (required).\n");
		fprintf (stderr, "-E sets file name for input estimated bathymetry model to polish.\n");
		fprintf (stderr, "Standard in is read for old-style imgbm pairs of (index,value) constraints.\n");
		fprintf (stderr, "img-style mercatorized grid coordinates only are supported.\n");
		fprintf (stderr, "You must specify -L72 or -L81 and -M1 or -M2 to indicate latitude range and pixel width of img file.\n"); 
		fprintf(stderr,"\tOPTIONS:\n");
		fprintf(stderr,"\t-C<c> set convergence limit to <c>, a positive even integer.  [%d]\n", CONVRG);
		fprintf(stderr,"\t-G<g> set initial grid mesh to every <g>'th point.  [Default finds GCF of nx and ny]\n");
		fprintf(stderr,"\t-I<i> set max number of iterations on each mesh to <i>.  [%d]\n", MAX_ITER);
		fprintf(stderr,"\t-T<t> set tension parameter to <t>.  [%4.2lf]\n", TENSION);
		exit(-1);
	}
	
	if (set_dimensions (&dims, &start_grid, latmax, minutes) ) {
		fprintf (stderr, "imgbm_polish:  WARNING.  Your initial grid mesh choice did not fit nx = %g and ny = %g.\n",
					(double)dims.nx, (double)dims.ny);
		fprintf (stderr, "\tIt was reset to %d, which is the greatest common factor of nx and ny.\n", start_grid);
	}
	
	if ( (hitmask = (unsigned char *) malloc (dims.nmask)) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot malloc for hitmask.\n");
		exit (EXIT_FAILURE);
	}
	
	if ( (seamask = (unsigned char *) malloc (dims.nmask)) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot malloc for seamask.\n");
		exit (EXIT_FAILURE);
	}
	
	if ( (ain = (short int *) malloc (2*dims.nxny)) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot malloc for array ain.\n");
		exit (EXIT_FAILURE);
	}
	
	if ( (a = (short int *) malloc (2*dims.nxmy)) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot malloc for array a.\n");
		exit (EXIT_FAILURE);
	}
	
	memset(a, 0, 2*dims.nxmy);	/* initialize the perturbation grid to zero */

	if (read_landsea_mask(bfile)) {
		fprintf(stderr,"imgbm_polish:  ERROR reading ocean/not-ocean 0/1 bitmask file %s\n", bfile);
		exit(EXIT_FAILURE);
	}
	

	f20 = -(1.0 - tension)/(20.0 - 16.0 * tension);
	f10 = (8.0 - 7.0*tension)/(20.0 - 16.0 * tension);
	f11 = 2.0 * f20;
	fprintf(stderr,"Factor 2,0:  %.5lg   Factor 1,1:  %.5lg   Factor 1,0:  %.5lg   Sum:  %.5lg\n", f20, f11, f10, 4.0*(f20+f11+f10));

	if (read_input_img(efile)) {
		fprintf(stderr,"imgbm_polish:  ERROR reading oldtopo file %s\n", efile);
		exit(-1);
	}

	if (read_constraints()) {
		fprintf(stderr,"imgbm_polish:  ERROR reading combined imgbm constraint data\n");
		exit(-1);
	}

	grid = start_grid * 5;
	last_grid = grid;
	n_tries = 0;
	do {
		grid = divide_grid (grid);
		if (n_tries) {
			interpolate_new_mesh(last_grid, grid);
		}
		else {
			initial_reset(grid);
		}
		iterate(grid, f10, f11, f20, max_iter, converge);
		last_grid = grid;
		n_tries++;
	} while (grid > 1);

	do_output ();

	exit(0);
}


int read_input_img(efile)
char	*efile;
{
	FILE *fp;
	int	i;
	short int emin = 0, emax = 0;

	if ( (fp = fopen(efile, "r")) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot open r input topo file %s\n", efile);
		return(-1);
	}

	if ( (fread((char *)ain, 2, dims.nxny, fp)) != dims.nxny) {
		fprintf (stderr, "imgbm_polish:  FATAL ERROR cannot read input topo file %s\n", efile);
		fclose (fp);
		return(-1);
	}

	fclose(fp);

	for (i = 0; i < dims.nxny; i++) {
		if (ain[i] < emin) emin = ain[i];
		if (ain[i] > emax) emax = ain[i];
	}

	fprintf(stderr,"imgbm_polish:  SUCCESS reading input img file.  Min = %hd  Max = %hd\n", emin, emax);

	return(0);
}

int	read_constraints()
{
	struct IMGPAIR dummy;
	int	n_read = 0;
	int	n_duplicate = 0;
	int	n_bad_index = 0;
	int	n_bad_value = 0;
	int	n_edit = 0;
	int	n_landsea_conflict = 0;
	int	ij, delta, max_delta = 0;
	int	ocean;

	/* First zero the constraint hit mask:  */
	memset(hitmask, 0, dims.nmask);

	while ( (fread((char *)&dummy, sizeof(struct IMGPAIR), 1, stdin)) == 1) {
	
		n_read++;

		if (dummy.index < 0 || dummy.index >= dims.nxny) {
			n_bad_index++;
			continue;
		}

		/* Might want to modify this step.  Used to eliminate extreme values:  */
		if (dummy.value > TALLEST || dummy.value < DEEPEST) {
			n_bad_value++;
			continue;
		}

		/* This step skips a point if we have already read a constraint for that point.
			The consequence is that when combining multiple constraint files together,
			they should be input in the order of precedence of the data.  */
		if ( (hitmask[dummy.index/8] & maskset[dummy.index%8]) != 0) {
			n_duplicate++;
			continue;
		}
		
		/* The next line sets ocean to 1 or 0, assuming seamask bits are ON when NOT ocean.  */
		ocean = ((seamask[dummy.index/8] & maskset[dummy.index%8]) == 0) ? 1 : 0;
		
		if (ocean && dummy.value >= 0) {
			n_landsea_conflict++;
			continue;
		}

		ij = dummy.index + dims.ijoff;
		delta = abs(dummy.value - ain[dummy.index]);
		if (delta > MAX_DELTA && dummy.value != -1 && dummy.value < 1000) {
			/* This edits large changes that are not at coastlines or high land topography */
			n_edit++;
			continue;
		}
		
		if (max_delta < delta) max_delta = delta;
		a[ij] = (short)(dummy.value - ain[dummy.index]);
		hitmask[dummy.index/8] |= maskset[dummy.index%8];
	}

	fprintf(stderr,"imgbm_polish:  STATUS REPORT.  Finished reading %d imgbm constraints.\n", n_read);
	fprintf(stderr,"The maximum difference between accepted constraints and the input img was %d\n", max_delta);
	fprintf(stderr,"%d inputs were skipped due to bad array indices.\n", n_bad_index);
	fprintf(stderr,"%d inputs were skipped because values exceeded %d or %d.\n", n_bad_value, DEEPEST, TALLEST);
	fprintf(stderr,"%d inputs were skipped because they would cause change > %d.\n", n_edit, MAX_DELTA);
	fprintf(stderr,"%d inputs were skipped because they violated the ocean mask.\n", n_landsea_conflict);
	fprintf(stderr,"%d inputs were skipped because they would constrain a pixel more than once.\n", n_duplicate);
	if (n_duplicate) fprintf(stderr,"WARNING:  Only the first valid constraint is used in case of conflict.  Run imgbm_combine or enter data in desired sequence.\n");

	return(0);
}


void iterate(grid, f10, f11, f20, max_iter, converge)
int	grid, max_iter, converge;
double	f10, f11, f20;
{
	int	i, j, ij, jgrid, iter_count, max_change, k;
	int	up2, up1, dn1, dn2, il2, il1, ir1, ir2;
	int	grid2, jgrid2, maxj, maxi;

	double	asum, adiff, sor;
	int	aadj, last_max_change;

	void	set_bc();

	iter_count = 0;
	jgrid = dims.nx * grid;
	jgrid2 = 2*jgrid;
	grid2 = 2*grid;
	sor = SOR;
	last_max_change = 5000;

	do {
		max_change = 0;
		iter_count++;

		set_bc(grid, jgrid, jgrid2);

		for (j = 0; j < dims.ny; j += grid) {
			if (j == 0) {
				up2 = -dims.ijoff;
				up1 = -dims.nx;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			else if (j == grid) {
				up2 = -jgrid - dims.nx;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			else if (j == dims.ny - grid) {
				up2 = -jgrid2;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid + dims.nx;
			}
			else {
				up2 = -jgrid2;
				up1 = -jgrid;
				dn1 = jgrid;
				dn2 = jgrid2;
			}
			for (i = 0; i < dims.nx; i+= grid) {
				ij = j * dims.nx + i;
				if (hitmask[ij/8] & maskset[ij%8]) continue;	/* This spot is known  */
				k = ij;
				ij += dims.ijoff;
				if (i == 0) {
					il2 = dims.nx - grid2;
					il1 = dims.nx - grid;
					ir1 = grid;
					ir2 = grid2;
				}
				else if (i == grid) {
					il2 = dims.nx - grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = grid2;
				}
				else if (i == dims.nx - grid2) {
					il2 = -grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = -i;
				}
				else if (i == dims.nx - grid) {
					il2 = -grid2;
					il1 = -grid;
					ir1 = -i;
					ir2 = -i + grid;
				}
				else {
					il2 = -grid2;
					il1 = -grid;
					ir1 = grid;
					ir2 = grid2;
				}

/*				if (debug) {
					if (ij + il2 < 0 || ij + il2 >= dims.nxmy) fprintf(stderr,"BUG on il2 at %d %d %d\n", grid, j, i);
					if (ij + ir2 < 0 || ij + ir2 >= dims.nxmy) fprintf(stderr,"BUG on ir2 at %d %d %d\n", grid, j, i);
					if (ij + up2 < 0 || ij + up2 >= dims.nxmy) fprintf(stderr,"BUG on up2 at %d %d %d\n", grid, j, i);
					if (ij + dn2 < 0 || ij + dn2 >= dims.nxmy) fprintf(stderr,"BUG on dn2 at %d %d %d\n", grid, j, i);

					if (ij + il1 < 0 || ij + il1 >= dims.nxmy) fprintf(stderr,"BUG on il1 at %d %d %d\n", grid, j, i);
					if (ij + ir1 < 0 || ij + ir1 >= dims.nxmy) fprintf(stderr,"BUG on ir1 at %d %d %d\n", grid, j, i);
					if (ij + up1 < 0 || ij + up1 >= dims.nxmy) fprintf(stderr,"BUG on up1 at %d %d %d\n", grid, j, i);
					if (ij + dn1 < 0 || ij + dn1 >= dims.nxmy) fprintf(stderr,"BUG on dn1 at %d %d %d\n", grid, j, i);

					if (ij+il1+up1 < 0 || ij+il1+up1 >= dims.nxmy) fprintf(stderr,"BUG on il1 up1 at %d %d %d\n", grid, j, i);
					if (ij+il1+dn1 < 0 || ij+il1+dn1 >= dims.nxmy) fprintf(stderr,"BUG on il1 dn1 at %d %d %d\n", grid, j, i);
					if (ij+ir1+up1 < 0 || ij+ir1+up1 >= dims.nxmy) fprintf(stderr,"BUG on ir1 up1 at %d %d %d\n", grid, j, i);
					if (ij+ir1+dn1 < 0 || ij+ir1+dn1 >= dims.nxmy) fprintf(stderr,"BUG on ir1 dn1 at %d %d %d\n", grid, j, i);
				}
*/


				asum = f20 * (a[ij + il2] + a[ij + ir2] + a[ij + up2] + a[ij + dn2])
					+ f10 * (a[ij + il1] + a[ij + ir1] + a[ij + up1] + a[ij + dn1])
					+ f11 * (a[ij+il1+up1]+a[ij+il1+dn1]+a[ij+ir1+up1]+a[ij+ir1+dn1]);
				if (rint(fabs(asum)) > 32000.0) fprintf(stderr,"BUG:  Range error at %d %d %d\n", grid, j, i);
				adiff = asum - a[ij];
				aadj = abs((int)rint(adiff * sor));
				/* Safety check here:  Don't change by more than 5000 m  */
				if (aadj > 5000) aadj = 5000;
				if (max_change < aadj) {
					max_change = aadj;
					maxj = j;
					maxi = i;
				}
				a[ij] = (adiff < 0.0) ? a[ij] - aadj : a[ij] + aadj;
				if (a[ij] + ain[k] >= 0 && (seamask[k/8] & maskset[k%8]) == 0) {
					/* This perturbation would pull us above sea level where we should stay below.
						Set the solution to MOST_SHOAL  */
					a[ij] = MOST_SHOAL - ain[k];
				}
			}
		}

		if (max_change > last_max_change) sor = 0.5*(1.0 + sor);

		last_max_change = max_change;

		fprintf(stderr,"%d\tD\t%d\t%d\t%d\t%d\n", grid, iter_count, max_change, maxj, maxi);

	} while (iter_count < max_iter && max_change > converge);

	return;
}


void	set_bc(grid, jgrid, jgrid2)
int	grid, jgrid, jgrid2;
{
	int	i, ij, il1, ir1;

	/* First set Top row d2a/dn2 = 0:  */
	for (i = 0; i < dims.nx; i+=grid) {
		a[dims.nx+i] = 2 * a[dims.ijoff+i] - a[dims.ijoff + jgrid + i];
	}

	/* Next set Top row d(Laplacian(a))/dn = 0:  */
	for (i = 0; i < dims.nx; i += grid) {
		il1 = (i == 0) ? dims.nx-grid : -grid;
		ir1 = (i == dims.nx-grid) ? grid-dims.nx : grid;
		a[i] = a[i+dims.ijoff+jgrid2] + a[i+dims.ijoff+jgrid+il1] + a[i+dims.ijoff+jgrid+ir1]
			-4*a[i+dims.ijoff+jgrid] +4*a[i+dims.nx] - a[i+dims.nx+il1] - a[i+dims.nx+ir1];
	}

	/* Next set Bottom row d2a/dn2 = 0:  */
	for (i = 0; i < dims.nx; i+=grid) {
		ij = dims.nxny + dims.ijoff + i;
		a[ij] = 2 * a[ij - jgrid] - a[ij - jgrid2];
	}

	/* Finally, set Bottom row d(Laplacian(a))/dn = 0:  */
	for (i = 0; i < dims.nx; i += grid) {
		il1 = (i == 0) ? dims.nx-grid : -grid;
		ir1 = (i == dims.nx-grid) ? grid-dims.nx : grid;
		ij = i + dims.nxny + dims.ijoff - jgrid;
		a[ij+dims.nx+jgrid] = a[ij-jgrid2] + a[ij-jgrid+il1] + a[ij-jgrid+ir1]
			-4*a[ij-jgrid] +4*a[ij+jgrid] - a[ij+jgrid+il1] - a[ij+jgrid+ir1];
	}

	return;
}	


void	initial_reset(grid)
int	grid;
{
	/* Visit every grid'th point and if it is not already constrained, set it equal to zero.  */
	int	i, j, ij;
	for (j = 0; j < dims.ny; j += grid) {
		for (i = 0; i < dims.nx; i+= grid) {
			ij = j * dims.nx + i;
			if (hitmask[ij/8] & maskset[ij%8]) continue;
			a[ij+dims.ijoff] = 0;
		}
	}
	return;
}

void	interpolate_new_mesh(old_grid, grid)
int	old_grid, grid;
{
	int	i, j, ij, ir, ii, jnx, j2, jj;
	double	dx;

	/* For each old row, interpolate new columns:  */
	for (j = 0; j < dims.ny; j += old_grid) {
		jnx = j * dims.nx;
		for (i = 0; i < dims.nx; i += old_grid) {
			ir = i + old_grid;
			if (ir == dims.nx) ir = 0;
			for (ii = i + grid; ii < i + old_grid; ii += grid) {
				ij = jnx + ii;
				if (hitmask[ij/8] & maskset[ij%8]) continue;
				dx = (double)(ii-i)/(double)old_grid;
				a[ij+dims.ijoff] = (int)rint(dx*a[jnx+ir+dims.ijoff] + (1.0 - dx)*a[jnx+i+dims.ijoff]);
			}
		}
	}

	/* Now go down all columns in new rows:  */
	for (j = 0; j < dims.ny; j+= old_grid) {
		j2 = j + old_grid;
		for (jj = j + grid; jj < j + old_grid; jj += grid) {
			dx = (double)(jj-j)/(double)old_grid;
			for (i = 0; i < dims.nx; i+= grid) {
				ij = jj * dims.nx + i;
				if (hitmask[ij/8] & maskset[ij%8]) continue;
				if (j2 == dims.ny) {
					a[ij+dims.ijoff] = 0;
				}
				else {
					a[ij+dims.ijoff] = (int)rint(dx*a[dims.ijoff+j2*dims.nx+i] + (1.0-dx)*a[dims.ijoff+j*dims.nx+i]);
				}
			}
		}
	}
	return;
}


void do_output (void)
{

	/* This routine does not [yet?] check whether the ocean mask is being violated.  */

	int	i, j, ij;

	for (j = 0; j < dims.ny; j++) {
		for (i = 0; i < dims.nx; i++) {
			ij = j * dims.nx + i;
			ain[ij] += a[ij + dims.ijoff];
			
			if (hitmask[ij/8] & maskset[ij%8]) {
				/* This point should be odd.  */
				if (ain[ij]%2 == 0) {
					if (ain[ij] >= 0)
						ain[ij]++;
					else
						ain[ij]--;
				}
			}
			else {
				/* This point should be even  */
				if (ain[ij]%2 != 0) {
					if (ain[ij] > 0)
						ain[ij]--;
					else
						ain[ij]++;
				}				
			}
		}
	}

	if ( (fwrite((char *)ain, 2, dims.nxny, stdout)) != dims.nxny) {
		fprintf(stderr,"Write failure sending solution to stdout.\n");
		exit(-1);
	}
	return;
}

int	set_dimensions (struct DIMS *dims, int *start_grid, int latmax, int minutes) {

	int	gcf_default, return_value = 0;
		
	if (latmax == 81) {
		dims->ny = 8640;
		gcf_default = 2160;
	}
	else {
		dims->ny = 6336;
		gcf_default = 144;
	}
	
	if (minutes == 1) {
		dims->ny *= 2;
		gcf_default *= 2;
		dims->nx = 21600;
	}
	else {
		dims->nx = 10800;
	}
	
	dims->nxny = dims->nx * dims->ny;
	dims->nxmy = dims->nx * (dims->ny + 4);	/* Two extra rows at top and bottom store boundary conditions  */
	dims->nmask = (dims->nxny)/8;	/* Number of bytes in each bitmask.  */
	dims->ijoff = dims->nx * 2;	/* Start of first data point, after skipping two top rows of BCs.  */
	
	if ((*start_grid) != 0) {
		/* User-supplied value must be tested.  */
		if (  (((dims->nx)%(*start_grid)) != 0) || (((dims->ny)%(*start_grid)) != 0) ) {
			*start_grid = gcf_default;
			return_value = -1;
		}
	}
	else {
		/* User did not supply a value.  Use default.  */
		*start_grid = gcf_default;
	}
	return (return_value);
}

int	divide_grid (int grid) {

	/* Grid is a positive integer greater than 1 when this routine is called.
		Grid = k, where every k'th point in the grid is being worked on.
		Thus k should be at most the greatest common factor of nx and ny.  For all
		supported coordinate sizes, k is of the form 2^a * 3^b * 5^c.  Thus
		this routine first divides grid by 5, if it is divisible by 5.  If
		not, we divide by 3, if it is divisible by 3.  if not, we divide by
		2.  This sequence (dividing by largest prime factor each time) was
		found to be best when "surface" was developed for GMT.
		This routine should not be called if grid is already equal to 1,
		or if it has some other peculiar number.  */
		
	if ( (grid%5) == 0) return (grid/5);
	if ( (grid%3) == 0) return (grid/3);
	if ( (grid%2) == 0) return (grid/2);
	return (grid);
}

int	read_landsea_mask (char *bfile) {

	FILE	*fp;
	
	if ( (fp = fopen (bfile, "r")) == NULL) {
		fprintf (stderr, "imgbm_polish:  FATAL_ERROR cannot open r %s\n", bfile);
		return (-1);
	}
	
	if ( (fread ((void *)seamask, (size_t)1, dims.nmask, fp)) != dims.nmask) {
		fprintf (stderr, "imgbm_polish:  FATAL_ERROR cannot read %lg bytes from %s\n", (double)dims.nmask, bfile);
		fclose (fp);
		return (-1);
	}
	fclose (fp);
	fprintf (stderr, "imgbm_polish:  SUCCESS reading landsea mask.\n");
	return (0);
}