/*----------------------------------------------------------------------*\
 | Program to check the integrity of .bna files produced by read_lgo	|
 |																		|
 | Usage: check_bna [-c] [-t deg] input_file							|
 |																		|
 | Options:	-c		close polygons if not closed						|
 |			-t deg	set tolerance to deg degrees; polygon vertices with	|
 |					external angles less than this value are removed.	|
 |																		|
 | Peter N. Schweitzer (USGS, Reston, VA 22092)							|
\*----------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#define MAXLEN 128

static double tolerance = 0.5;

static int memicmp (const char *s, const char *t, int n) {
	int r = 0;
	while (n-- > 0 && (r = toupper (*s++) - toupper (*t++)) == 0);
	return (r);
	}

static int purify (double *x, double *y, int n);

main (int argc, char *argv[]) {
	int i,j,k,m,n;
	int np;
	FILE *in;
	char input_file[MAXLEN];
	char string[MAXLEN],*s;
	char primary_id[17],label[64],*p;
	double *x,*y;
	int close = 0;

	if (argc > 1) {
		for (i=0; i < argc; i++) {
			if (memcmp (argv[i],"-c",2) == 0) close = 1;
			else
				if (memcmp (argv[i],"-t",2) == 0) {
					i++;
					tolerance = atof (argv[i]);
					}
				else
					strcpy (input_file,argv[i]);
			}
		}
	else {
		fprintf (stderr,"Usage: %s [-c] [-t deg] infile\n",argv[0]);
		exit (1);
		}

	if (in = fopen (input_file,"r")) {
		k = 0;
		while (fgets (string,MAXLEN,in)) {
			k++;
			if (s = strrchr (string,'\n')) *s = 0;

			p = primary_id;
			for (j=0, s=string; *s && *s != '\t'; j++, s++, p++) if (j < 16) *p = *s;
			*p = 0;
			if (*s == 0) {
				fprintf (stderr,"Error (line %d): no point count for \"%s\"\n",k,primary_id);
				exit (1);
				}

			while (*s && isspace (*s)) s++;
			p = label;
			for (j=0; *s && *s != '\t'; j++, p++, s++) if (j < 64) *p = *s;
			*p = 0;
			if (*s == 0) {
				fprintf (stderr,"Error (line %d): no point count for \"%s\"\n",k,primary_id);
				exit (1);
				}
			while (*s && isspace (*s)) s++;

			if (np = atoi (s)) {
				fprintf (stderr,"%s\t\"%s\"\t%d\t",primary_id,label,np);
				if (!(x = (double *) malloc ((np+1) * sizeof (double)))) {
					fprintf (stderr,"Error: could not allocate point array\n");
					exit (1);
					}
				if (!(y = (double *) malloc ((np+1) * sizeof (double)))) {
					fprintf (stderr,"Error: could not allocate point array\n");
					exit (1);
					}
				for (j=0; j < np; j++) {
					fgets (string,MAXLEN,in);
					k++;
					if (sscanf (string,"%lf%lf",&x[j],&y[j]) != 2) {
						fprintf (stderr,"Error (line %d): two doubles expected in \"%s\"\n",k,string);
						exit (1);
						}
					}
				if (close)
					if (x[0] != x[np-1] || y[0] != y[np-1]) {
						fprintf (stderr,"(closing)\t");
						x[np] = x[0];
						y[np] = y[0];
						np++;
						}
				n = purify (x,y,np);
				fprintf (stderr,"%d\n",np - n);
				printf ("%s\t%s\t%d\n",primary_id,label,n);
				for (j=0; j < n; j++) printf ("%.6lf\t%.6lf\n",x[j],y[j]);
				free (x);
				free (y);
				}
			else {
				fprintf (stderr,"%s\t\"%s\"\t%d\n",primary_id,label,np);
				}
			}
		fclose (in);
		}
	else {
		fprintf (stderr,"Error: could not open input file %s\n",input_file);
		exit (1);
		}
	}

static int purify (double *x, double *y, int n) {
	int i,j,k,m;
	int *q;
	double ax,ay,bx,by,a,b,*t;
	double small,d2r;

	d2r = 3.141592653589793/180.0;

	small = cos(d2r * tolerance);

	if (!(q = (int *) malloc (n * sizeof (int)))) {
		fprintf (stderr,"Error: could not allocate q array\n");
		exit (1);
		}
	if (!(t = (double *) malloc (n * sizeof (double)))) {
		fprintf (stderr,"Error: could not allocate t array\n");
		exit (1);
		}

	for (i=0; i < n; i++) q[i] = 1;

	/*------------------------------------------------------------------*\
	 | Eliminate duplicate points										|
	\*------------------------------------------------------------------*/

	i = 0;
	for (j=1; j < n; j++)
		if (x[j] == x[i] && y[j] == y[i]) q[j] = 0;
		else i = j;

	i = 0;
	for (j=0; j < n; j++)
		if (q[j]) {
			x[i] = x[j];
			y[i] = y[j];
			q[i] = 1;
			i++;
			}
	n = i;

	/*------------------------------------------------------------------*\
	 | Eliminate intermediate points on lines							|
	\*------------------------------------------------------------------*/

	for (j=1; j < n-1; j++) {
		i = j-1;
		k = j+1;
		ax = x[j] - x[i];
		ay = y[j] - y[i];
		bx = x[k] - x[j];
		by = y[k] - y[j];
		a = ax*ax + ay*ay;
		b = bx*bx + by*by;
		t[j] = fabs(ax*bx + ay*by)/sqrt(a*b);
		if (t[j] > small) q[j] = 0;		/* remember cos(small) is near 1 */
		}

	i = 0;
	for (j=0; j < n; j++)
		if (q[j]) {
			x[i] = x[j];
			y[i] = y[j];
			q[i] = 1;
			i++;
			}
	n = i;

	/*------------------------------------------------------------------*\
	 | Return the number of points retained.							|
	\*------------------------------------------------------------------*/

	return (n);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/