/* cavestats.c */
/* Print statistics on a cave */
/* Copyright (C) 2004 David Loeffler
 * Based on 3dtopos.c by Olly Betts
 *
 * 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; either version 2 of the License, or
 * (at your option) 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

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

#include "cmdline.h"
#include "filelist.h"
#include "filename.h"
#include "img.h"
#include "message.h"
#include "osalloc.h"

int cmp_station(const void *a, const void *b);

static const struct option long_opts[] = {
	/* const char *name; int has_arg (0 no_argument, 1 required_*, 2 optional_*); int *flag; int val; */
	{"survey", required_argument, 0, 's'},
	{"help", no_argument, 0, HLP_HELP},
	{"version", no_argument, 0, HLP_VERSION},
	{"compact", no_argument, 0, 'c'},
	{0, 0, 0, 0}
};

#define short_opts "s:c"

static struct help_msg help[] = {
/*				<-- */
	{HLP_ENCODELONG(0),			"only load the sub-survey with this prefix"},
	{0, 0}
};

typedef struct {
	const char *name;
	img_point pt;
} station;

static station basept;

static int separator;

int
cmp_station(const void *a, const void *b)
{
	img_point pta, ptb;
	double crossprod;
	pta = ((const station *)a)->pt;
	ptb = ((const station *)b)->pt;

	if(pta.x == basept.pt.x && pta.y == basept.pt.y) return -1;
	if(ptb.x == basept.pt.x && ptb.y == basept.pt.y) return 1;
	/* These two lines mean that the basepoint should compare as being before everything.
	Otherwise somewhat unpredictable behaviour occurs! */
	crossprod = (pta.x - basept.pt.x)*(ptb.y - basept.pt.y) - (pta.y - basept.pt.y)*(ptb.x - basept.pt.x);
	if(crossprod > 0) return 1;
	if(crossprod < 0) return -1;
	else return 0;
}

static station *stns;
static OSSIZE_T c_stns = 0;
static station *hull;
static OSSIZE_T c_hull = 0;

int
main(int argc, char **argv)
{
	char *fnm;
	img *pimg;
	int result;
	OSSIZE_T c_labels = 0;
	OSSIZE_T c_totlabel = 0;
	char *p, *p_end;
	const char *survey = NULL;
	bool compact_output = 0;
	double totlen = 0, planlen = 0;
	img_point lastpt;
	int num_legs = 0;
	
	msg_init(argv);

	cmdline_set_syntax_message("3D_FILE", NULL); /* TRANSLATE */
	cmdline_init(argc, argv, short_opts, long_opts, NULL, help, 1, 1);
	while (1) {
		int opt = cmdline_getopt();
		if (opt == EOF) break;
		if (opt == 's') survey = optarg;
		if (opt == 'c') compact_output = 1;
	}

	fnm = argv[optind++];

	pimg = img_open_survey(fnm, survey);
	if (!pimg) fatalerror(img_error(), fnm);
	separator = pimg->separator;

	/* Output headings line */
	if (!compact_output) printf("Statistics for %s (prefix %s)\n", fnm, survey);
	
	do {
		img_point pt;
		result = img_read_item(pimg, &pt);
		switch (result) {
			case img_LINE:
					if(!pimg->flags)
					{
						planlen += sqrt((pt.x - lastpt.x)*(pt.x - lastpt.x) + (pt.y - lastpt.y)*(pt.y - lastpt.y));
						totlen += sqrt((pt.x - lastpt.x)*(pt.x - lastpt.x) + (pt.y - lastpt.y)*(pt.y - lastpt.y) + (pt.z - lastpt.z)*(pt.z - lastpt.z));
					}
					num_legs++;
			case img_MOVE:
					lastpt = pt;
				break;
			case img_LABEL:
				if(pimg->flags & img_SFLAG_UNDERGROUND)
				{
					c_labels++;
					c_totlabel += strlen(pimg->label) + 1;
				}
				break;
			case img_BAD:
				img_close(pimg);
				fatalerror(img_error(), fnm);
		}
	} while (result != img_STOP);

	/* + 1 to allow for reading coordinates of legs after the last label */
	stns = osmalloc(ossizeof(station) * (c_labels + 1));
	p = osmalloc(c_totlabel);
	p_end = p + c_totlabel;

	img_rewind(pimg);

	do {
	  result = img_read_item(pimg, &(stns[c_stns].pt));
	  switch (result) {
		case img_LINE:
		case img_MOVE:
	 break;
		case img_LABEL:
	if (!(pimg->flags & img_SFLAG_UNDERGROUND)) break; 
	if (c_stns < c_labels) {
		OSSIZE_T len = strlen(pimg->label) + 1;
		if (p + len <= p_end) {
			memcpy(p, pimg->label, len);
			stns[c_stns++].name = p;
			p += len;
			break;
		}
	 }
	 /* FALLTHRU */
		case img_BAD:
	 img_close(pimg);
	 fatalerror(/*Bad 3d image file `%s'*/106, fnm);
	  }
	} while (result != img_STOP);

	if (c_stns != c_labels || p != p_end) {
	  img_close(pimg);
	  fatalerror(/*Bad 3d image file `%s'*/106, fnm);
	}

	img_close(pimg);
	if (c_stns == 0)
	{
		printf("No survey data!\n");
		return EXIT_SUCCESS;
	}
	station wmost, nmost, emost, smost, umost, lmost;
	wmost = emost = nmost = smost = umost = lmost = stns[0];
	{
		OSSIZE_T i;
		for(i = 0; i < c_stns; i++)
		{
			if(stns[i].pt.x < wmost.pt.x) wmost = stns[i];
			if(stns[i].pt.x > emost.pt.x) emost = stns[i];
			if(stns[i].pt.y > nmost.pt.y) nmost = stns[i];
			if(stns[i].pt.y < smost.pt.y) smost = stns[i];
			if(stns[i].pt.z > umost.pt.z) umost = stns[i];
			if(stns[i].pt.z < lmost.pt.z) lmost = stns[i]; 
		}
	}
	
	basept = wmost; /* any extreme point will do */
	
	qsort(stns, c_stns, sizeof(station), cmp_station);

	/* Now, having sorted the stations, any duplicate stations caused by *equates will appear consecutive. */
	
	{
			OSSIZE_T i;
			int real_stn_count = 1;
			for(i = 1; i < c_stns; i++)
					if(stns[i].pt.x != stns[i-1].pt.x || stns[i].pt.y != stns[i-1].pt.y || stns[i].pt.z != stns[i-1].pt.z) 
							real_stn_count++;
			if (!compact_output) 
			{
				printf((real_stn_count == 1 ? msg(172) : msg(173)), real_stn_count);
				printf((num_legs == 1 ? msg(174) : msg(175)), num_legs);
				putnl();
				//printf((num_legs - real_stn_count == 0 ? msg(138) : msg(139)), num_legs - real_stn_count + 1);
				// Can't do this, as we don't know the number of connected components
			}
	}
	

	if (!compact_output)
	{
		printf("Total length of survey legs = %.2fm\n", totlen);

		printf(msg(133), planlen);
		putnl();
	
		printf(msg(135), (umost.pt.z - lmost.pt.z));
		printf(umost.name);
		printf(msg(136), umost.pt.z);
		printf(lmost.name);
		printf(msg(137), lmost.pt.z);
		putnl();
	
		printf(msg(148), (nmost.pt.y - smost.pt.y));
		printf(nmost.name);
		printf(msg(136), nmost.pt.y);
		printf(smost.name);
		printf(msg(137), smost.pt.y);
		putnl();
	
		printf(msg(149), (emost.pt.x - wmost.pt.x));
		printf(emost.name);
		printf(msg(136), emost.pt.x);
		printf(wmost.name);
		printf(msg(137), wmost.pt.x);
		putnl();
	}
	if(c_stns == 1)
	{
		if (!compact_output) printf("Horizontal extent = 0.00m (from %s to %s)\n", stns[0].name, stns[0].name);
		else printf("0.00\t0.00\t0.00\n");
		return EXIT_SUCCESS;
	}
	
	if(c_stns == 2)
	{
		if (!compact_output) printf("Horizontal extent = %.2fm (from %s to %s)\n", planlen, stns[0].name, stns[1].name);
		else printf("%.2f\t%.2f\t%.2f", totlen, umost.pt.z-lmost.pt.z, planlen);
		return EXIT_SUCCESS;
	}

	// Noddy equivalent of a stack...
	c_hull = 3;
	hull = osmalloc(ossizeof(station) * c_stns);
	hull[0] = stns[0];
	hull[1] = stns[1];
	hull[2] = stns[2];
	OSSIZE_T i = 3;
	while(i < c_stns)
	{
		double sa;
		sa = (hull[c_hull - 1].pt.x - stns[i].pt.x)*(hull[c_hull - 2].pt.y - stns[i].pt.y)
				- (hull[c_hull - 2].pt.x - stns[i].pt.x)*(hull[c_hull - 1].pt.y - stns[i].pt.y);
		if(sa > 0)
		{
			hull[c_hull++] = stns[i++];
		}
		else
		{
			c_hull--;
		}	
	}
	{
		OSSIZE_T j;
		station s1, s2;
		double d, e, r;
		s1 = s2 = hull[0];
		r = 0;
		j = 0;
		d = 0;
		for(i = 0; i < c_hull; i++)
		{
				d = sqrt((hull[j % c_hull].pt.x - hull[i].pt.x)*(hull[j % c_hull].pt.x - hull[i].pt.x) + (hull[j % c_hull].pt.y - hull[i].pt.y)*(hull[j % c_hull].pt.y - hull[i].pt.y));
				while((e = sqrt((hull[(j+1) % c_hull].pt.x - hull[i].pt.x)*(hull[(j+1) % c_hull].pt.x - hull[i].pt.x) + (hull[(j+1) % c_hull].pt.y - hull[i].pt.y)*(hull[(j+1) % c_hull].pt.y - hull[i].pt.y))) >= d)
				{
						j++; d = e;
				}
				if(d > r)
				{
						r = d;
						s1 = hull[i];
						s2 = hull[j % c_hull];
				}
		}
		if (!compact_output) printf("Horizontal extent = %.2fm (from %s to %s)\n", r, s1.name, s2.name);
		else printf("%.2f\t%.2f\t%.2f\n", totlen, umost.pt.z - lmost.pt.z, r);
	}
	return EXIT_SUCCESS;
}