/*
 * File: trian.c
 *
 * (c) Peter Kleiweg
 *     2006
 *
 * This 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,
 * or (at your option) any later version.
 *
 */

#define my_VERSION "0.02"

#define __NO_MATH_INLINES

#ifdef __WIN32__
#  define my_PATH_SEP '\\'
#else
#  define my_PATH_SEP '/'
#endif

#include <ctype.h>
#include <errno.h>
#include <float.h>
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>

typedef enum { FALSE = 0, TRUE = 1} BOOL;

typedef struct {
    int
        i;
    char
        *s;
    float
        x,
        y;
    BOOL
        defined;
} NAME;

typedef struct point_ {
    int
        i;
    double
        x,
        y;
    struct point_
        *next;
} POINT;

NAME
    *names = NULL;

POINT
    **point;

#define BUFSIZE 2048

BOOL
    asDiff = FALSE;

char
    **arg_v,
    *configfile,
    *outfile = NULL,
    *labelfile = NULL,
    *coordfile = NULL,
    buffer [BUFSIZE + 1],
    buf    [BUFSIZE + 1],
    *programname,
    *no_mem_buffer,
    out_of_memory [] = "Out of memory!";

int
    arg_c,
    n_names = 0,
    max_names = 0,
    high_names = 0,
    n_lines = 0,
    max_lines = 0,
    volgende = 1,
    top,
    inputline;

double
    margin = 1.0e-20,
    aspect = 1.0,
    xp,
    yp;

time_t
    tp;

FILE
    *fp_out;

void
    crop (int, int),
    get_programname (char const *argv0),
    process_args (void),
    errit (char const *format, ...),
    warn (char const *format, ...),
    syntax (void),
    *s_malloc (size_t size),
    *s_realloc (void *block, size_t size);

char
    *get_arg (void),
    *s_strdup (char const *s);

char const
    *unquote (char *s);

int
    icmp (void const *p1, void const *p2),
    scmp (void const *p1, void const *p2),
    lblcmp (void const *key, void const *p);

BOOL
    getline (FILE *fp);

FILE
    *r_fopen (char const *filename);

BOOL EQUAL (double d1, double d2)
{
    return (d1 - d2 <= margin && d2 - d1 <= margin) ? TRUE : FALSE;
}

BOOL ZERO (double d)
{
    return EQUAL (d, 0.0);
}

BOOL RANGE (double d, double d1, double d2, BOOL incl)
{
    if (incl) {
        if ((d >= d1 && d <= d2) || (d >= d2 && d <= d1))
            return TRUE;
        if (EQUAL (d, d1) || EQUAL (d, d2))
            return TRUE;
        return FALSE;
    }

    if (EQUAL (d, d1) || EQUAL (d, d2))
        return FALSE;
    if ((d > d1 && d < d2) || (d > d2 && d < d1))
        return TRUE;
    return FALSE;
}

int main (int argc, char *argv [])
{
    int
        i,
	j,
        j1,
        j2,
        n,
	**dm;
    double
        x,
        y,
        dx,
        dy,
        xmin,
        xmax,
        ymin,
        ymax;
    FILE
        *fp;
    POINT
        *ptmp;
    NAME
        *n_ptr;

    no_mem_buffer = (char *) malloc (1024);

    get_programname (argv [0]);

    arg_c = argc;
    arg_v = argv;
    process_args ();
    
    if (arg_c != 2)
        syntax ();    

    coordfile = arg_v [1];

    if (outfile) {
        fp_out = fopen (outfile, "w");
        if (! fp_out)
            errit ("Creating file \"%s\": %s", outfile, strerror (errno));
    } else
        fp_out = stdout;

    if (labelfile) {
        fp = r_fopen (labelfile);
        while (getline (fp)) {
            if (sscanf (buffer, "%i%n", &i, &n) < 1)
                errit ("Missing index number in file \"%s\", line %i", labelfile, inputline);
            if (i < 1)
                errit ("Illegal value in file \"%s\", line %i", labelfile, inputline);
            while (buffer [n] && isspace ((unsigned char) buffer [n]))
                n++;    
            if (! buffer [n])
                errit ("Missing name in file \"%s\", line %i", labelfile, inputline);
            if (n_names == max_names) {
                max_names += 256;
                names = (NAME *) s_realloc (names, max_names * sizeof (NAME));
            }
            names [n_names].i = i;
            names [n_names].s = s_strdup (unquote (buffer + n));
            names [n_names++].defined = FALSE;
        }
        fclose (fp);
        qsort (names, n_names, sizeof (NAME), icmp);
        for (i = 0; i < n_names; i++) {
            if (names [i].i < i + 1)
                errit ("Duplicate labels found for number %i in file \"%s\"", names [i].i, labelfile);
            if (names [i].i > i + 1)
                errit ("Missing label for number %i in file \"%s\"", i + 1, labelfile);
        }
        qsort (names, n_names, sizeof (NAME), scmp);	
    }


    fp = r_fopen (coordfile);
    while (getline (fp)) {
        if (sscanf (buffer, "%lf %lf %lf %lf%n", &x, &y, &dx, &dy, &n) < 4)
            errit ("Missing coordinates in file \"%s\", line %i", arg_v [1], inputline);
        while (buffer [n] && isspace ((unsigned char) buffer [n]))
            n++;
        if (! buffer [n])
            errit ("Missing name in file \"%s\", line %i", arg_v [1], inputline);
        if (labelfile) {
            n_ptr = bsearch (unquote (buffer + n), names, n_names, sizeof (NAME), lblcmp);
            if (! n_ptr)
                continue;
        } else {
            if (max_names == n_names) {
                max_names += 256;
                names = (NAME *) s_realloc (names, max_names * sizeof (NAME));
            }
            names [n_names].s = s_strdup (unquote (buffer + n));
            n_ptr = names + n_names++;
        }
        n_ptr->x = x;
        n_ptr->y = y;
        n_ptr->defined = TRUE;
    }
    fclose (fp);

    if (labelfile) {
        qsort (names, n_names, sizeof (NAME), icmp);
        for (i = 0; i < n_names; i++)
            if (! names [i].defined)
                errit ("Missing coordinates for \"%s\"", names [i].s);
    }

    xmin = ymin = FLT_MAX;
    xmax = ymax = -FLT_MAX;
    for (i = 0; i < n_names; i++) {
	names [i].y *= aspect;
	if (names [i].x < xmin)
	    xmin = names [i].x;
	if (names [i].x > xmax)
	    xmax = names [i].x;
	if (names [i].y < ymin)
	    ymin = names [i].y;
	if (names [i].y > ymax)
	    ymax = names [i].y;
    }
    point = (POINT **) s_malloc (n_names * sizeof (POINT *));
    for (i = 0; i < n_names; i++) {
	point [i] = (POINT *) s_malloc (sizeof (POINT));
	point [i]->x = xmin - (xmax - xmin);
	point [i]->y = ymin - (ymax - ymin);
	point [i]->next = NULL;
	point [i]->i = -1;

	ptmp = (POINT *) s_malloc (sizeof (POINT));
	ptmp->i = -1;
	ptmp->x = xmax + (xmax - xmin);
	ptmp->y = ymin - (ymax - ymin);
	ptmp->next = point [i];
	point [i] = ptmp;

	ptmp = (POINT *) s_malloc (sizeof (POINT));
	ptmp->i = -1;
	ptmp->x = xmax + (xmax - xmin);
	ptmp->y = ymax + (ymax - ymin);
	ptmp->next = point [i];
	point [i] = ptmp;

	ptmp = (POINT *) s_malloc (sizeof (POINT));
	ptmp->i = -1;
	ptmp->x = xmin - (xmax - xmin);
	ptmp->y = ymax + (ymax - ymin);
	ptmp->next = point [i];
	point [i] = ptmp;
    }
    for (i = 0; i < n_names; i++) {
	if (((n_names - i) % 10) == 0)
	    fprintf (stderr, "\r%i ", n_names - i);
	fflush (stderr);
	for (j1 = i + 1, j2 = i - 1; j1 < n_names || j2 >= 0; j1++, j2--) {
	    if (j1 < n_names) {
		if (EQUAL (names [i].x, names [j1].x) && EQUAL (names [i].y, names [j1].y))
		    warn ("Identical coordinates for \"%s\" and \"%s\"", names [i].s, names [j1].s);
		else
		    crop (i, j1);
	    }
	    if (j2 >= 0)
		if (! (EQUAL (names [i].x, names [j2].x) && EQUAL (names [i].y, names [j2].y)))
		    crop (i, j2);
	}
    }
    fprintf (stderr, "\r    \r");
    for (i = 0; i < n_names; i++) {
	names [i].y /= aspect;
	for (ptmp = point [i]; ptmp != NULL; ptmp = ptmp->next)
	    ptmp->y /= aspect;
    }

    fprintf (fp_out, "%i\n", n_names);

    for (i = 0; i < n_names; i++)
	fprintf (fp_out, "%s\n", names [i].s);

    if (asDiff) {
	dm = (int **) s_malloc (n_names * sizeof (int *));
	for (i = 0; i < n_names; i++) {
	    dm [i] = (int *) s_malloc (n_names * sizeof (int));
	    for (j = 0; j < n_names; j++)
		dm [i][j] = 0;
	}
	for (i = 0; i < n_names; i++)
	    for (ptmp = point [i]; ptmp != NULL; ptmp = ptmp->next)
		if (ptmp->i > i)
		    dm [i][ptmp->i] = dm [ptmp->i][i] = 1;
	for (i = 1; i < n_names; i++)
	    for (j = 0; j < i; j++)
		fprintf (fp_out, "%i\n", dm [i][j]);
    } else {
	for (i = 0; i < n_names; i++)
	    for (ptmp = point [i]; ptmp != NULL; ptmp = ptmp->next)
		if (ptmp->i > i)
		    fprintf (fp_out, "%i %i\n", i + 1, ptmp->i + 1);
    }

    if (outfile)
        fclose (fp_out);

    return 0;
}

/*
 * (x1, y1) and (x2, y2) are endpoints of short line
 * (x, y) is point on long line
 * (rx, ry) is direction of long line
 * result in (xp, yp)
 */
BOOL docross (double x1, double y1, double x2, double y2, double x, double y, double rx, double ry, BOOL incl)
{
    double
        a,
        aa,
        b,
        bb,
        rrx,
        rry;

    if (ZERO (rx) && ZERO (ry))
        return FALSE;

    /* direction of short line */
    rrx = x2 - x1;
    rry = y2 - y1;

    if (ZERO (rrx) && ZERO (rry))
        return FALSE;

    /* long line vertical */
    if (ZERO (rx)) {

        /* both lines vertical */
        if (ZERO (rrx))
            return FALSE;

        xp = x;
        a = (x - x1) / rrx;
        yp = y1 + a * rry;
        return RANGE (a, 0.0, 1.0, incl);
    }

    /* short line vertical */
    if (ZERO (rrx)) {
        xp = x1;
        a = (x1 - x) / rx;
        yp = y + a * ry;
        return RANGE (yp, y1, y2, incl);
    }

    a  = ry  / rx;
    aa = rry / rrx;
    if (EQUAL (a, aa))    /* parallel lines */
        return FALSE;
    b  = y  - a  * x;
    bb = y1 - aa * x1;
    xp = (bb - b) / (a - aa);
    yp = a * xp + b;
    return RANGE (xp, x1, x2, incl);
}

void crop (int i, int j)
{
    double
        x,
        y,
        rx,
        ry,
        xp1,
        yp1,
        xp2,
        yp2;
    POINT
        *p,
        *p1,
        *p2,
        *p3,
        *p4;

    x = (names [i].x + names [j].x) / 2.0;
    y = (names [i].y + names [j].y) / 2.0;
    rx = names [j].y - names [i].y;
    ry = names [i].x - names [j].x;

    for (p1 = point [i]; p1->next != NULL; p1 = p1->next) {
        if (p1 == p1->next)
            errit ("p1 == p1->next: THIS IS A BUG");
        p2 = p1->next;
        if (docross (p1->x, p1->y, p2->x, p2->y, x, y, rx, ry, TRUE))
            break;
    }
    if (! p1->next)
         return;
    xp1 = xp;
    yp1 = yp;

    for (p3 = p2; p3 != NULL; p3 = p3->next) {
        if (p3 == p3->next)
            errit ("p3 == p3->next: THIS IS A BUG");
        p4 = p3->next;
        if (! p4)
            p4 = point [i];
        if (docross (p3->x, p3->y, p4->x, p4->y, x, y, rx, ry, TRUE))
            if (! (EQUAL (xp1, xp) && EQUAL (yp1, yp)))
                break;
    }
    if (! p3)
        return;
    xp2 = xp;
    yp2 = yp;

    /* p1 on line */
    if (EQUAL (p1->x, xp1) && EQUAL (p1->y, yp1)) {
        if (docross (p2->x, p2->y, names [i].x, names [i].y, xp1, yp1, xp2 - xp1, yp2 - yp1, FALSE)) {
            /* p2 at outside */
	    p1->i = j;
            if (EQUAL (xp2, p4->x) && EQUAL (yp2, p4->y))
                p1->next = p4;
            else {
                p1->next = p3;
                p3->x = xp2;
                p3->y = yp2;
            }
        } else {
            /* p2 at inside */
            if (EQUAL (xp2, p3->x) && EQUAL (yp2, p3->y)) {
		p3->i = j;
                p3->next = NULL;
            } else if (p3->next) {
		p4->i = j;
                p4->x = xp2;
                p4->y = yp2;
                p4->next = NULL;
            }
        }
        return;
    }

    /* p2 on line */
    if (EQUAL (p2->x, xp1) && EQUAL (p2->y, yp1)) {
        if (p2 == p3)
            return;
        if (docross (p1->x, p1->y, names [i].x, names [i].y, xp1, yp1, xp2 - xp1, yp2 - yp1, FALSE)) {
            /* p1 at outside */
            point [i] = p2;
            if (! p3->next)
                p3->next = p4 = (POINT *) s_malloc (sizeof (POINT));
	    p4->i = j;
            p4->x = xp2;
            p4->y = yp2;
            p4->next = NULL;
        } else {
            /* p1 at inside */
            if (EQUAL (xp2, p3->x) && EQUAL (yp2, p3->y))
                return;
	    p2->i = j;
            if (EQUAL (xp2, p4->x) && EQUAL (yp2, p4->y))
                p2->next = p3->next;
            else {
                p2->next = p3;
                p3->x = xp2;
                p3->y = yp2;
            }
        }
        return;
    }

    if (docross (p1->x, p1->y, names [i].x, names [i].y, xp1, yp1, xp2 - xp1, yp2 - yp1, FALSE)) {
        /* p1 at outside */
        point [i] = p1;
        p1->x = xp1;
        p1->y = yp1;
        if (! p3->next)
            p4 = p3->next = (POINT *) s_malloc (sizeof (POINT));
	p4->i = j;
        p4->x = xp2;
        p4->y = yp2;
        p4->next = NULL;
    } else {
        /* p1 at inside */
        p2->x = xp1;
        p2->y = yp1;
        if (EQUAL (p4->x, xp2) && EQUAL (p4->y, yp2))
            p2->next = p4;
        else {
            if (p2 == p3) {
                p = p3->next;
                p3 = (POINT *) s_malloc (sizeof (POINT));
                p3->next = p;
		p3->i = p2->i;
            }
            p2->next = p3;
            p3->x = xp2;
            p3->y = yp2;
        }
	p2->i = j;
    }
}

int icmp (void const *p1, void const *p2)
{
    return ((NAME *)p1)->i - ((NAME *)p2)->i;
}

int scmp (void const *p1, void const *p2)
{
    return strcmp (((NAME *)p1)->s, ((NAME *)p2)->s);
}

int lblcmp (void const *key, void const *p)
{
    return strcmp ((char *) key, ((NAME *)p)->s);
}


char const *unquote (char *s)
{
    int
        i,
        j;

    if (s [0] != '"')
        return s;
    j = 0;
    for (i = 1; s [i]; i++) {
        if (s [i] == '"')
            break;
        if (s [i] == '\\') {
            if (s [i + 1])
                s [j++] = s [++i];
        } else
            s [j++] = s [i];
    }
    s [j] = '\0';
    return s;
}

FILE *r_fopen (char const *filename)
{
    FILE
        *fp;

    inputline = 0;
    fp = fopen (filename, "r");
    if (! fp)
        errit ("Opening file \"%s\": %s", filename, strerror (errno));
    return fp;
}

BOOL getline (FILE *fp)
/* Lees een regel in
 * Plaats in buffer
 * Negeer lege regels en regels die beginnen met #
 */
{
    int
        i;

    for (;;) {
        if (fgets (buffer, BUFSIZE, fp) == NULL)
            return FALSE;

        inputline++;
        i = strlen (buffer);
        while (i && isspace ((unsigned char) buffer [i - 1]))
            buffer [--i] = '\0';
        i = 0;
        while (buffer [i] && isspace ((unsigned char) buffer [i]))
            i++;
        if (buffer [i] == '#')
            continue;
        if (buffer [i]) {
            memmove (buffer, buffer + i, strlen (buffer) + 1);
            return TRUE;
        }
    }
}

void process_args ()
{
    while (arg_c > 1 && arg_v [1][0] == '-') {
        switch (arg_v [1][1]) {
	    case 'a':
		aspect = atof (get_arg ());
		if (aspect < 1.0)
		    errit ("Illegal value for 'aspect' in file \"%s\", line %i", configfile, inputline);
		break;
	    case 'd':
		asDiff = TRUE;
		break;
            case 'l':
                labelfile = get_arg ();
                break;
            case 'o':
                outfile = get_arg ();
                break;
            default:
                errit ("Illegal option '%s'", arg_v [1]);
        }
        arg_c--;
        arg_v++;
    }
}

char *get_arg ()
{
    if (arg_v [1][2])
        return arg_v [1] + 2;

    if (arg_c == 2)
        errit ("Missing argument for '%s'", arg_v [1]);

    arg_v++;
    arg_c--;
    return arg_v [1];
}

void errit (char const *format, ...)
{
    va_list
        list;

    fprintf (stderr, "\nError %s: ", programname);

    va_start (list, format);
    vfprintf (stderr, format, list);

    fprintf (stderr, "\n\n");

    exit (1);
}

void warn (char const *format, ...)
{
    va_list
        list;

    fprintf (stderr, "\nWarning %s: ", programname);

    va_start (list, format);
    vfprintf (stderr, format, list);

    fprintf (stderr, "\n\n");
}

void get_programname (char const *argv0)
{
#ifdef __MSDOS__
    char
        name [MAXFILE];
    fnsplit (argv0, NULL, NULL, name, NULL);
    programname = strdup (name);
#else
    char
        *p;
    p = strrchr (argv0, my_PATH_SEP);
    if (p)
        programname = strdup (p + 1);
    else
        programname = strdup (argv0);
#endif
}

void *s_malloc (size_t size)
{
    void
        *p;

    p = malloc (size);
    if (! p) {
        free (no_mem_buffer);
        errit (out_of_memory);
    }
    return p;
}

void *s_realloc (void *block, size_t size)
{
    void
        *p;

    p = realloc (block, size);
    if (! p) {
        free (no_mem_buffer);
        errit (out_of_memory);
    }
    return p;
}

char *s_strdup (char const *s)
{
    char
        *s1;

    if (s) {
        s1 = (char *) s_malloc (strlen (s) + 1);
        strcpy (s1, s);
    } else {
        s1 = (char *) s_malloc (1);
        s1 [0] = '\0';
    }
    return s1;
}

void syntax ()
{
    fprintf (
        stderr,
        "\n"
        "Triangulation, Version " my_VERSION "\n"
        "(c) P. Kleiweg 2006\n"
        "\n"
        "Usage: %s [-a float] [-d] [-l filename] [-o filename] coordinate_file\n"
	"\n"
	"\t-a : aspect ratio, (>= 1)\n"
	"\t-d : output in format of a difference matrix file\n"
	"\t-o : output file\n"
        "\n",
        programname
    );
    exit (1);
}