/*--------------------------------------------------------------------
 *	$Id: triangulate.c,v 1.3.4.3 2002/02/27 17:58:55 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith
 *	See COPYING 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 General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	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 General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * triangulate reads one or more files (or GMT_stdin) with x,y[,whatever] and
 * outputs the indeces of the vertices of the optimal Delaunay triangulation
 * using the method by Watson, D. F., ACORD: Automatic contouring of raw data,
 * Computers & Geosciences, 8, 97-101, 1982.  Optionally, the output may take
 * the form of (1) a multi-segment file with the vertice coordinates needed to
 * draw the triangles, or (2) a grdfile based on gridding the plane estimates
 * PS. Instead of Watson's method you may choose to link with the triangulate
 * routine written by Jonathan Shewchuck.  See the file TRIANGLE.HOWTO for
 * details.  That function is far faster than Watson's method.
 *
 * Author:      Paul Wessel
 * Date:        1-JAN-1993
 * Version:     2.0
 * Revised:	13-AUG-1998 for GMT 3.1
 * Revised:	02-MAR-1999 for GMT 3.2
 * Revised:	29-MAR-2000 for GMT 3.3.5
 * Modified:	10 Jul 2000 by PW to add -L
 * Version:	3.4.1
 *
 */
 
#include "gmt.h"

struct EDGE {
	int begin, end;
};
	
main (int argc, char **argv)
{
	int i, j, ij, ij1, ij2, ij3, np, k, fno, n = 0, n_alloc = 0, n_files = 0, n_edge;
	int n_expected_fields, n_args, *link, n_fields, compare_edge(const void *p1, const void *p2);
	int one_or_zero, i_min, i_max, j_min, j_max, p;
	
	BOOLEAN error = FALSE, map_them = FALSE, inc_set = FALSE, do_ddx = FALSE, do_ddy = FALSE;
	BOOLEAN nofile = TRUE, done = FALSE, do_grid = FALSE, set_empty = FALSE;
	
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, zj, zk, zl, zlj, zkj, *in;
	double *xx, *yy, *zz, vx[4], vy[4], xinc2, yinc2, idx, idy, xp, yp, a, b, c, f;
	double xkj, xlj, ykj, ylj, out[2];
	
	float *grd, empty = 0.0;
	
	char line[BUFSIZ], *outfile = CNULL;
	
	FILE *fp = NULL;
	
	struct GRD_HEADER header;

	struct EDGE *edge;
	
	argc = GMT_begin (argc, argv);

	GMT_grd_init (&header, argc, argv, FALSE);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
        
				/* Common parameters */
        
				case 'H':
				case 'J':
				case 'R':
				case 'V':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
        
				/* Supplemental parameters */
        
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'D':
					do_ddx = (argv[i][2] == 'x' || argv[i][2] == 'X');
					do_ddy = (argv[i][2] == 'y' || argv[i][2] == 'Y');
					break;
				case 'E':
					empty = (argv[i][2] == 'N' || argv[i][2] == 'n') ? GMT_f_NaN : (float)atof (&argv[i][2]);
					set_empty = TRUE;
					break;
				case 'F':
					header.node_offset = TRUE;
					break;
				case 'G':
					outfile = &argv[i][2];
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
					inc_set = TRUE;
					break;
				case 'L':
					GMT_geographic_in = TRUE;
					break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "triangulate %s - Optimal (Delaunay) triangulation of xyz-data\n\n", GMT_VERSION);
		fprintf (stderr, "usage: triangulate <infiles> [-Dx|y] [-E<empty>] [-F] [-G<grdfile> [-H[<nrec>]] [-I<dx>[m|c][/<dy>[m|c]]] [-J<parameters>]\n");
		fprintf (stderr, "\t[-L] [-M[<flag>]] [-R<west/east/south/north>] [-V] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
                
		if (GMT_quick) exit (EXIT_FAILURE);
                
		fprintf (stderr, "\tinfiles (in ASCII) has 2 or more columns.  If no file(s) is given, standard input is read.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-Dx or -Dy takes derivative in that direction (only with -G) [Default is z value].\n");
		fprintf (stderr, "\t-E value to use for empty nodes [Default is NaN].\n");
		fprintf (stderr, "\t-F Force pixel registration (only with -G) [Default is gridline registration]\n");
		fprintf (stderr, "\t-G Grid data. Give name of output gridfile and specify -R -I\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I sets the grid spacing for the grid.  Append m for minutes, c for seconds\n");
		GMT_explain_option ('J');   
		fprintf (stderr, "\t-L means that x is longitude, i.e. assumed to be periodic in 360\n");
		fprintf (stderr, "\t-M output triangle edges as multiple segments separated by a record\n");
		fprintf (stderr, "\t   whose first character is <flag> ['>']\n");
		fprintf (stderr, "\t   [Default is to output the indeces of vertices for each Delaunay triangle]\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		fprintf (stderr, "\t-bo writes binary integer index table; ignored if -M is selected [Default is ASCII i/o]\n");
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (GMT_io.binary[0] && gmtdefs.io_header) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
        if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < 2) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 2 columns\n", GMT_program);
		error++;
	}
	
	do_grid = (outfile && inc_set && header.x_inc > 0.0 && header.y_inc > 0.0);
	
	if (do_grid && !project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Must specify -R, -I, -G for gridding\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */

	if (GMT_io.binary[0] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	if (do_grid) {
		header.x_min = west;	header.x_max = east;
		header.y_min = south;	header.y_max = north;
		GMT_grd_RI_verify (&header, 1);
	}

	if (west != east && project_info.projection > 0) { /* Gave -R -J */
		map_them = TRUE;
		GMT_map_setup (west, east, south, north);
	}

	/* Now we are ready to take on some input values */
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2 + do_grid;
	
	n_alloc = GMT_CHUNK;
	xx = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	yy = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	if (do_grid) zz = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
		
	n = 0;
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over input files, if any */
		if (!nofile && argv[fno][0] == '-') continue;
		
		if (nofile) {	/* Just read standard input */
			fp = GMT_stdin;
			done = TRUE;
#ifdef SET_IO_MODE
			GMT_setmode (0);
#endif
		}
		else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]);
			continue;
		}

		if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Reading file %s\n", GMT_program, argv[fno]);
		
		if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp);

		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */
					
			if (GMT_io.status & GMT_IO_MISMATCH) {
				fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n);
				exit (EXIT_FAILURE);
			}

			xx[n] = in[0];
			yy[n] = in[1];
			if (do_grid) zz[n] = in[2];
			n++;
			
			if (n == n_alloc) {	/* Get more memory */
				n_alloc += GMT_CHUNK;
				xx = (double *) GMT_memory ((void *)xx, (size_t)n_alloc, sizeof (double), GMT_program);
				yy = (double *) GMT_memory ((void *)yy, (size_t)n_alloc, sizeof (double), GMT_program);
				if (do_grid) zz = (double *) GMT_memory ((void *)zz, (size_t)n_alloc, sizeof (double), GMT_program);
			}
		}
		
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	xx = (double *) GMT_memory ((void *)xx, (size_t)n, sizeof (double), GMT_program);
	yy = (double *) GMT_memory ((void *)yy, (size_t)n, sizeof (double), GMT_program);
	if (do_grid) zz = (double *) GMT_memory ((void *)zz, (size_t)n, sizeof (double), GMT_program);

	if (map_them) {	/* Must make parallel arrays for projected x/y */

		double *xxp, *yyp;

		xxp = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
		yyp = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
		for (i = 0; i < n; i++) GMT_geo_to_xy (xx[i], yy[i], &xxp[i], &yyp[i]);

		if (gmtdefs.verbose) fprintf (stderr, "%s: Do Delaunay optimal triangulation on projected coordinates\n", GMT_program);

		np = GMT_delaunay (xxp, yyp, n, &link);

		GMT_free ((void *)xxp);
		GMT_free ((void *)yyp);
	}
	else {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Do Delaunay optimal triangulation on given coordinates\n", GMT_program);

		np = GMT_delaunay (xx, yy, n, &link);
	}

	if (gmtdefs.verbose) fprintf (stderr, "%s: %d Delaunay triangles found\n", GMT_program, np);

	if (do_grid) {
	
		if (header.node_offset) {
			one_or_zero = 0;
			xinc2 = 0.5 * header.x_inc;
			yinc2 = 0.5 * header.y_inc;
		}
		else {
			one_or_zero = 1;
			xinc2 = yinc2 = 0.0;
		}
		idx = 1.0 / header.x_inc;
		idy = 1.0 / header.y_inc;
		header.nx = irint ( (header.x_max - header.x_min) * idx) + one_or_zero;
		header.ny = irint ( (header.y_max - header.y_min) * idy) + one_or_zero;
		grd = (float *) GMT_memory (VNULL, (size_t)(header.nx * header.ny), sizeof (float), GMT_program);
		if (!set_empty) empty = GMT_f_NaN;
		for (i = 0; i < header.nx * header.ny; i++) grd[i] = empty;	/* initialize grid */
		
		for (k = ij = 0; k < np; k++) {
		
			/* Find equation for the plane as z = ax + by + c */
			
			vx[0] = vx[3] = xx[link[ij]];	vy[0] = vy[3] = yy[link[ij]];	zj = zz[link[ij++]];
			vx[1] = xx[link[ij]];	vy[1] = yy[link[ij]];	zk = zz[link[ij++]];
			vx[2] = xx[link[ij]];	vy[2] = yy[link[ij]];	zl = zz[link[ij++]];
			
			xkj = vx[1] - vx[0];	ykj = vy[1] - vy[0];
			zkj = zk - zj;	xlj = vx[2] - vx[0];
			ylj = vy[2] - vy[0];	zlj = zl - zj;
			
			f = 1.0 / (xkj * ylj - ykj * xlj);
			a = -f * (ykj * zlj - zkj * ylj);
			b = -f * (zkj * xlj - xkj * zlj);
			c = -a * vx[1] - b * vy[1] + zk;

			/* Compute grid indices the current triangle may cover, assuming all triangles are
			   in the -R region (header.x_min/x_max etc.)  Always, i_min <= i_max, j_min <= j_max.
			 */

			i_min = (int)floor ((MIN (MIN (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
			i_max = (int)ceil ((MAX (MAX (vx[0], vx[1]), vx[2]) - header.x_min + xinc2) * idx);
			j_min = (int)floor ((header.y_max - MAX (MAX (vy[0], vy[1]), vy[2]) + yinc2) * idy);
			j_max = (int)ceil ((header.y_max - MIN (MIN (vy[0], vy[1]), vy[2]) + yinc2) * idy);
			
			/* Adjustments for triangles outside -R region. */
			/* Triangle to the left or right. */
			if ((i_max < 0) || (i_min >= header.nx)) continue;
			/* Triangle Above or below */
			if ((j_max < 0) || (j_min >= header.ny)) continue;

			/* Triangle covers boundary, left or right. */
			if (i_min < 0) i_min = 0;       if (i_max >= header.nx) i_max = header.nx - 1;
			/* Triangle covers boundary, top or bottom. */
			if (j_min < 0) j_min = 0;       if (j_max >= header.ny) j_max = header.ny - 1;

			for (j = j_min; j <= j_max; j++) {
				yp = header.y_max - j * header.y_inc - yinc2;
				p = j * header.nx + i_min;
				for (i = i_min; i <= i_max; i++, p++) {
					xp = header.x_min + i * header.x_inc + xinc2;
					
					if (!GMT_non_zero_winding (xp, yp, vx, vy, 4)) continue;	/* Outside */
					
					if (do_ddx)
						grd[p] = (float)a;
					else if (do_ddy)
						grd[p] = (float)b;
					else
						grd[p] = (float)(a * xp + b * yp + c);
				}
			}
		}
		if (GMT_write_grd (outfile, &header, grd, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
			exit (EXIT_FAILURE);
		}
					
	}
	if (GMT_io.multi_segments) {	/* Must find unique edges to output only once */
		n_edge = 3 * np;
		edge = (struct EDGE *) GMT_memory (VNULL, (size_t)n_edge, sizeof (struct EDGE), GMT_program);
		for (i = ij1 = 0, ij2 = 1, ij3 = 2; i < np; i++, ij1 += 3, ij2 += 3, ij3 += 3) {
			edge[ij1].begin = link[ij1];	edge[ij1].end = link[ij2];
			edge[ij2].begin = link[ij2];	edge[ij2].end = link[ij3];
			edge[ij3].begin = link[ij1];	edge[ij3].end = link[ij3];
		}
		for (i = 0; i < n_edge; i++) if (edge[i].begin > edge[i].end) i_swap (edge[i].begin, edge[i].end);
		
		qsort ((void *)edge, (size_t)n_edge, sizeof (struct EDGE), compare_edge);
		for (i = 1, j = 0; i < n_edge; i++) {
			if (edge[i].begin != edge[j].begin || edge[i].end != edge[j].end) j++;
			edge[j] = edge[i];
		}
		n_edge = j + 1;
	
		if (gmtdefs.verbose) fprintf (stderr, "%s: %d unique triangle edges\n", GMT_program, n_edge);
		
		for (i = 0; i < n_edge; i++) {
			fprintf (GMT_stdout, "%c Edge %d-%d\n", GMT_io.EOF_flag, edge[i].begin, edge[i].end);
			out[0] = xx[edge[i].begin];	out[1] = yy[edge[i].begin];
			GMT_output (GMT_stdout, 2, out);
			out[0] = xx[edge[i].end];	out[1] = yy[edge[i].end];
			GMT_output (GMT_stdout, 2, out);
		}
		GMT_free ((void *)edge);
	}
	else if (GMT_io.binary[1])
		fwrite ((void *)link, sizeof (int), (size_t)(3*np), GMT_stdout);
	else
		for (i = ij = 0; i < np; i++, ij += 3) fprintf (GMT_stdout, "%d\t%d\t%d\n", link[ij], link[ij+1], link[ij+2]);
	
	GMT_free ((void *) xx);
	GMT_free ((void *) yy);
	GMT_free ((void *) link);
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program);
	
	GMT_end (argc, argv);
}

int compare_edge (const void *p1, const void *p2)
{
	struct EDGE *a, *b;

	a = (struct EDGE *)p1;
	b = (struct EDGE *)p2;
	if (a->begin < b->begin)
		return (-1);
	else if (a->begin > b->begin)
		return (1);
	else {
		if (a->end < b->end)
			return (-1);
		else if (a->end > b->end)
			return (1);
		else
			return (0);
	}
}