/* 3DVto2D.c  -  convert 3DV data to 2-D list format */

/* Oscar Garcia <garciao@mof.govt.nz>, May 1992 */

/* uses getopt */
int getopt(int argc, char *argv[], char *options);
extern  int optind;
extern char *optarg;

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

#define USAGE "usage: 3DVto2D [/tx] [/nx] [/dx] [input_file [output_file]]\n\
\tIf file names are omitted, the standard i/o streams are used.\n\
\tOptions (with defaults):\n\
\t\t/t0   -  Latitude, degrees\n\
\t\t/n0   -  Longitude, degrees\n\
\t\t/d-3  -  Viewing distance.  Zero gives orthogonal projection.\n\
\t\t             Positive is in user units to object center.\n\
\t\t             Negative is relative to object size.\n\
\t\t             Default -3 is as in 3dv."

#define ERROR(msg) {fputs(msg,stderr),exit(1);}
#define PERROR(msg) {perror(msg),exit(1);}

#define BIG 3e33
#define EPS 1e-9
#define RADIANS 3.141592/180
#define DDEFAULT -3

void main(int argc, char* argv[])
{
	int i, j, n;
	double scale, slon, clon, slat, clat, t[3][3],
		lat = 0, lon = 0, dist = DDEFAULT;
	char options[] = "t:n:d:";
	FILE *infile = stdin, *outfile = stdout;
	typedef double POINT[3];
	POINT min, max, centre, xx, *x;


	/* get options */
	while ((n = getopt(argc, argv, options)) != EOF)
		if (n == '?')
			ERROR(USAGE)
		else
			switch (n)
			{	case 't':	lat = atof(optarg); break;
				case 'n':	lon = atof(optarg); break;
				case 'd':	dist = atof(optarg); break;
			}

	/* open files */
	if (optind < argc - 2)
		ERROR(USAGE);	/* too many */
	if (optind < argc)
	{	infile = fopen(argv[optind], "rt");
		if (infile == NULL)
			ERROR("Can't open input file");
		if (++optind < argc)
		{	outfile = fopen(argv[optind], "wt");
			if (outfile == NULL)
				ERROR("Can't open output file");
		}
	}

	/* read points */
	i = fscanf(infile, "%d", &n);
	if (i != 1 || ferror(infile))
		PERROR("Bad input file");
#ifdef __MSDOS__
	if ((unsigned)n * sizeof(POINT) > 0xfff0U)
		ERROR("64k memory block limit exceeded");
#endif
	x = malloc((unsigned)n * sizeof(POINT));
	if (x == NULL)
		ERROR("Out of memory");
	for (j = 0; j < 3; j++)
	{	min[j] = BIG;
		max[j] = -BIG;
	}
	for (i = 0; i < n; i++)
	{	fscanf(infile, "%lf %lf %lf", &x[i][0], &x[i][1], &x[i][2]);
		for (j = 0; j < 3; j++)
		{	if (x[i][j] < min[j])
				min[j] = x[i][j];
			else if (x[i][j] > max[j])
				max[j] = x[i][j];
		}
	}
	if (ferror(infile))
		PERROR("Error reading input");

	/* get centre and distance */
	for (j = 0, scale = 0.0; j < 3; j++)
	{	centre[j] = (min[j] + max[j]) / 2;
		scale = scale + (max[j] - min[j]) * (max[j] - min[j]);
	}
	if (dist < -EPS)	/* perspective, relative dist; convert to absolute */
		dist *= sqrt(scale * (1 + 0.25 / (dist * dist))) / -0.99;

	/* transformation matrix */
	slon = sin(lon * RADIANS);
	clon = cos(lon * RADIANS);
	slat = sin(lat * RADIANS);
	clat = cos(lat * RADIANS);
	t[0][0] = clon;
	t[0][1] = - slon;
	t[1][0] = clat * slon;
	t[1][1] = clat * clon;
	t[1][2] = slat;
	t[2][0] = - slat * slon;
	t[2][1] = - slat * clon;
	t[2][2] = clat;

	/* rotate and project along y */
	min[0] = min[2] = BIG;	/* 2-d coords will be saved in 0 and 2 */
	max[0] = max[2] = -BIG;
	for (i = 0; i < n; i++)
	{	for (j = 0; j < 3; j++)
			xx[j] = x[i][j] - centre[j];
		x[i][0] = t[0][0] * xx[0] + t[0][1] * xx[1];
		x[i][1] = t[1][0] * xx[0] + t[1][1] * xx[1] + t[1][2] * xx[2];
		x[i][2] = t[2][0] * xx[0] + t[2][1] * xx[1] + t[2][2] * xx[2];
		if (dist > EPS)		/* perspective */
		{	x[i][0] /= 1 - x[i][1] / dist;
			x[i][2] /= 1 - x[i][1] / dist;
		}	/* x[i][1] not used again here, but might need it some day */
		if (x[i][0] < min[0])
			min[0] = x[i][0];
		else if (x[i][0] > max[0])
			max[0] = x[i][0];
		if (x[i][2] < min[2])
			min[2] = x[i][2];
		else if (x[i][2] > max[2])
			max[2] = x[i][2];
	}

	/* factor for scaling to [-1, 1] */
	scale = max[0];
	if (-min[0] > scale)
		scale = -min[0];
	if (max[2] > scale)
		scale = max[2];
	if (-min[2] > scale)
		scale = -min[2];
	scale = 1 / scale;

	/* output info */
	fscanf(infile, "%d", &n);
	fprintf(outfile, "# n=%d scale=%g x_range=[%g:%g] y_range=[%g:%g]\n",
		n, 1/scale, scale*min[0], scale*max[0], scale*min[2], scale*max[2]);
	fprintf(outfile, "# latitude=%g longitude=%g distance=%g\n",
		lat, lon, dist);

	/* output data */
	while (n-- > 0)
	{	fscanf(infile, "%d %d", &i, &j);
		if (j == 0)
			putc('\n', outfile);		/* blank line (move) */
		fprintf(outfile, "%g %g %d\n",
			scale * x[i-1][0], scale * x[i-1][2], j);
	}
	if (ferror(infile))
		PERROR("Error reading input");
	if (feof(infile))
		ERROR("Unexpected end of file on input");
	fclose(infile);
	if (ferror(outfile))
		PERROR("Error during output");
	fclose(outfile);

}