/*************************************************************************
 * Copyright (c) 2011 AT&T Intellectual Property 
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v1.0
 * which accompanies this distribution, and is available at
 * https://www.eclipse.org/legal/epl-v10.html
 *
 * Contributors: Details at https://graphviz.org
 *************************************************************************/

#include "config.h"

#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cgraph/cgraph.h>     /* for agerr() and friends */
#include <neatogen/delaunay.h>
#include <util/alloc.h>
#include <util/sort.h>

#ifdef HAVE_GTS
#include <gts.h>

static int triangle_is_hole(void *triangle, void *ignored) {
    GtsTriangle *t = triangle;
    (void)ignored;

    GtsEdge *e1, *e2, *e3;
    GtsVertex *v1, *v2, *v3;

    gts_triangle_vertices_edges(t, NULL, &v1, &v2, &v3, &e1, &e2, &e3);

    if ((GTS_IS_CONSTRAINT(e1) && GTS_SEGMENT(e1)->v1 != v1) ||
	(GTS_IS_CONSTRAINT(e2) && GTS_SEGMENT(e2)->v1 != v2) ||
	(GTS_IS_CONSTRAINT(e3) && GTS_SEGMENT(e3)->v1 != v3))
	return TRUE;

    return FALSE;
}

static unsigned delaunay_remove_holes(GtsSurface *surface) {
    return gts_surface_foreach_face_remove(surface, triangle_is_hole, NULL);
}

/* Derived classes for vertices and faces so we can assign integer ids
 * to make it easier to identify them. In particular, this allows the
 * segments and faces to refer to vertices using the order in which
 * they were passed in. 
 */
typedef struct {
    GtsVertex v;
    int idx;
} GVertex;

/// convert a `GtsVertex` pointer to a pointer to its derived class
static GVertex *downcast(GtsVertex *base) {
  return (GVertex *)((char *)base - offsetof(GVertex, v));
}

/// convert a `GVertex` pointer to a pointer to its base class
static GtsVertex *upcast(GVertex *derived) {
  return &derived->v;
}

typedef struct {
    GtsVertexClass parent_class;
} GVertexClass;

static GtsVertexClass *g_vertex_class(void) {
    static GVertexClass *klass = NULL;

    if (klass == NULL) {
	GtsObjectClassInfo vertex_info = {
	    .name = "GVertex",
	    .object_size = sizeof(GVertex),
	    .class_size = sizeof(GVertexClass),
	};
	klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_vertex_class()),
				     &vertex_info);
    }

    return &klass->parent_class;
}

typedef struct {
    GtsFace v;
    int idx;
} GFace;

typedef struct {
    GtsFaceClass parent_class;
} GFaceClass;

static GtsFaceClass *g_face_class(void) {
    static GFaceClass *klass = NULL;

    if (klass == NULL) {
	GtsObjectClassInfo face_info = {
	    .name = "GFace",
	    .object_size = sizeof(GFace),
	    .class_size = sizeof(GFaceClass),
	};
	klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_face_class()),
				     &face_info);
    }

    return &klass->parent_class;
}

/// destroy each edge using v, then destroy v
static void
destroy (GtsVertex* v)
{
    GSList * i;

    i = v->segments;
    while (i) {
	GSList * next = i->next;
	gts_object_destroy (i->data);
	i = next;
    }
    g_assert (v->segments == NULL);
    gts_object_destroy(GTS_OBJECT(v));
}

/* Main entry point to using GTS for triangulation.
 * Input is npt points with x and y coordinates stored either separately
 * in x[] and y[] (sepArr != 0) or consecutively in x[] (sepArr == 0).
 * Optionally, the input can include nsegs line segments, whose endpoint
 * indices are supplied in segs[2*i] and segs[2*i+1] yielding a constrained
 * triangulation.
 *
 * The return value is the corresponding gts surface, which can be queries for
 * the triangles and line segments composing the triangulation.
 */
static GtsSurface*
tri(double *x, double *y, int npt, int *segs, int nsegs, int sepArr)
{
    int i;
    GtsSurface *surface;
    GVertex **vertices = gv_calloc(npt, sizeof(GVertex *));
    GtsEdge **edges = gv_calloc(nsegs, sizeof(GtsEdge *));
    GSList *list = NULL;
    GtsVertex *v1, *v2, *v3;
    GtsTriangle *t;
    GtsVertexClass *vcl = g_vertex_class();
    GtsEdgeClass *ecl = GTS_EDGE_CLASS (gts_constraint_class ());

    if (sepArr) {
	for (i = 0; i < npt; i++) {
	    GVertex *p = downcast(gts_vertex_new(vcl, x[i], y[i], 0));
	    p->idx = i;
	    vertices[i] = p;
	}
    }
    else {
	for (i = 0; i < npt; i++) {
	    GVertex *p = downcast(gts_vertex_new(vcl, x[2*i], x[2*i+1], 0));
	    p->idx = i;
	    vertices[i] = p;
	}
    }

    /* N.B. Edges need to be created here, presumably before the
     * the vertices are added to the face. In particular, they cannot
     * be added created and added vi gts_delaunay_add_constraint() below.
     */
    for (i = 0; i < nsegs; i++) {
	edges[i] = gts_edge_new(ecl,
		 upcast(vertices[segs[2 * i]]),
		 upcast(vertices[segs[2 * i + 1]]));
    }

    for (i = 0; i < npt; i++)
	list = g_slist_prepend(list, vertices[i]);
    t = gts_triangle_enclosing(gts_triangle_class(), list, 100.);
    g_slist_free(list);

    gts_triangle_vertices(t, &v1, &v2, &v3);

    surface = gts_surface_new(gts_surface_class(), g_face_class(),
                              gts_edge_class(), gts_vertex_class());
    gts_surface_add_face(surface, gts_face_new(gts_face_class(),
					       t->e1, t->e2, t->e3));

    for (i = 0; i < npt; i++) {
	GtsVertex *v4 = upcast(vertices[i]);
	GtsVertex *v = gts_delaunay_add_vertex(surface, v4, NULL);

	/* if v != NULL, it is a previously added pt with the same
	 * coordinates as v4, in which case we replace v4 with v
	 */
	if (v && v4 != v) {
	    gts_vertex_replace(v4, v);
	}
    }

    for (i = 0; i < nsegs; i++) {
	gts_delaunay_add_constraint(surface,GTS_CONSTRAINT(edges[i]));
    }

    /* destroy enclosing triangle */
    gts_allow_floating_vertices = TRUE;
    gts_allow_floating_edges = TRUE;
    destroy(v1);
    destroy(v2);
    destroy(v3);
    gts_allow_floating_edges = FALSE;
    gts_allow_floating_vertices = FALSE;

    if (nsegs)
	delaunay_remove_holes(surface);

    free (edges);
    free(vertices);
    return surface;
}

typedef struct {
    int n;
    v_data *delaunay;
} estats;
    
static int cnt_edge(void *edge, void *stats) {
    GtsSegment *e = edge;
    estats *sp = stats;

    sp->n++;
    if (sp->delaunay) {
	sp->delaunay[downcast(e->v1)->idx].nedges++;
	sp->delaunay[downcast(e->v2)->idx].nedges++;
    }

    return 0;
}

static void
edgeStats (GtsSurface* s, estats* sp)
{
    gts_surface_foreach_edge(s, cnt_edge, sp);
}

static int add_edge(void *edge, void *data) {
    GtsSegment *e = edge;
    v_data *delaunay = data;

    int source = downcast(e->v1)->idx;
    int dest = downcast(e->v2)->idx;

    delaunay[source].edges[delaunay[source].nedges++] = dest;
    delaunay[dest].edges[delaunay[dest].nedges++] = source;

    return 0;
}

static v_data *delaunay_triangulation(double *x, double *y, int n) {
    GtsSurface* s = tri(x, y, n, NULL, 0, 1);
    int i, nedges;
    estats stats;

    if (!s) return NULL;

    v_data *delaunay = gv_calloc(n, sizeof(v_data));

    for (i = 0; i < n; i++) {
	delaunay[i].ewgts = NULL;
	delaunay[i].nedges = 1;
    }

    stats.n = 0;
    stats.delaunay = delaunay;
    edgeStats (s, &stats);
    nedges = stats.n;
    int *edges = gv_calloc(2 * nedges + n, sizeof(int));

    for (i = 0; i < n; i++) {
	delaunay[i].edges = edges;
	edges += delaunay[i].nedges;
	delaunay[i].edges[0] = i;
	delaunay[i].nedges = 1;
    }
    gts_surface_foreach_edge(s, add_edge, delaunay);

    gts_object_destroy (GTS_OBJECT (s));

    return delaunay;
}

typedef struct {
    int n;
    int* edges;
} estate;

static int addEdge(void *edge, void *state) {
    GtsSegment *e = edge;
    estate *es = state;

    int source = downcast(e->v1)->idx;
    int dest = downcast(e->v2)->idx;

    es->edges[2 * es->n] = source;
    es->edges[2 * es->n + 1] = dest;
    es->n += 1;

    return 0;
}

static int vcmp(const void *x, const void *y, void *values) {
    const int *a = x;
    const int *b = y;
    const double *_vals = values;
    double va = _vals[*a];
    double vb = _vals[*b];

    if (va < vb) return -1;
    if (va > vb) return 1;
    return 0;
}

/* Given n points whose coordinates are in the x[] and y[]
 * arrays, compute a Delaunay triangulation of the points.
 * The number of edges in the triangulation is returned in pnedges.
 * The return value itself is an array e of 2*(*pnedges) integers,
 * with edge i having points whose indices are e[2*i] and e[2*i+1].
 *
 * If the points are collinear, GTS fails with 0 edges.
 * In this case, we sort the points by x coordinates (or y coordinates
 * if the points form a vertical line). We then return a "triangulation"
 * consisting of the n-1 pairs of adjacent points.
 */
int *delaunay_tri(double *x, double *y, int n, int* pnedges)
{
    GtsSurface* s = tri(x, y, n, NULL, 0, 1);
    int* edges;

    if (!s) return NULL;

    estats stats = {0};
    edgeStats (s, &stats);
    int nedges = *pnedges = stats.n;

    if (nedges) {
	edges = gv_calloc(2 * nedges, sizeof(int));
	gts_surface_foreach_edge(s, addEdge, &(estate){.edges = edges});
    }
    else {
	int* vs = gv_calloc(n, sizeof(int));

	*pnedges = nedges = n-1;
	int *ip = edges = gv_calloc(2 * nedges, sizeof(int));

	for (int i = 0; i < n; i++)
	    vs[i] = i;

	gv_sort(vs, n, sizeof(int), vcmp, x[0] == x[1] /* vertical line? */ ? y : x);

	int tl = vs[0];
	for (int i = 1; i < n; i++) {
	    const int hd = vs[i];
	    *ip++ = tl;
	    *ip++ = hd;
	    tl = hd;
	}

	free (vs);
    }

    gts_object_destroy (GTS_OBJECT (s));

    return edges;
}

static int cntFace(void *face, void *data) {
    GFace *fp = face;
    int *ip = data;

    fp->idx = *ip;
    *ip += 1;

    return 0;
}

typedef struct {
    GtsSurface* s;
    int* faces;
    int* neigh;
} fstate;

typedef struct {
    int nneigh;
    int* neigh;
} ninfo;

static int addNeighbor(void *face, void *ni) {
    GFace *f = face;
    ninfo *es = ni;

    es->neigh[es->nneigh] = f->idx;
    es->nneigh++;

    return 0;
}

static int addFace(void *face, void *state) {
    GFace *f = face;
    fstate *es = state;

    int i, myid = f->idx;
    int* ip = es->faces + 3*myid;
    int* neigh = es->neigh + 3*myid;
    ninfo ni;
    GtsVertex *v1, *v2, *v3;

    gts_triangle_vertices (&f->v.triangle, &v1, &v2, &v3);
    *ip++ = downcast(v1)->idx;
    *ip++ = downcast(v2)->idx;
    *ip++ = downcast(v3)->idx;

    ni.nneigh = 0;
    ni.neigh = neigh;
    gts_face_foreach_neighbor((GtsFace*)f, 0, addNeighbor, &ni);
    for (i = ni.nneigh; i < 3; i++)
	neigh[i] = -1;

    return 0;
}

static int addTri(void *face, void *state) {
    GFace *f = face;
    fstate *es = state;

    int myid = f->idx;
    int* ip = es->faces + 3*myid;
    GtsVertex *v1, *v2, *v3;

    gts_triangle_vertices (&f->v.triangle, &v1, &v2, &v3);
    *ip++ = downcast(v1)->idx;
    *ip++ = downcast(v2)->idx;
    *ip++ = downcast(v3)->idx;

    return 0;
}

/* Given n points whose coordinates are in x[] and y[], and nsegs line
 * segments whose end point indices are given in segs, return a surface
 * corresponding the constrained Delaunay triangulation.
 * The surface records the line segments, the triangles, and the neighboring
 * triangles.
 */
surface_t* 
mkSurface (double *x, double *y, int n, int* segs, int nsegs)
{
    GtsSurface* s = tri(x, y, n, segs, nsegs, 1);
    estats stats;
    estate state;
    fstate statf;
    int nfaces = 0;

    if (!s) return NULL;

    surface_t *sf = gv_alloc(sizeof(surface_t));
    stats.n = 0;
    stats.delaunay = NULL;
    edgeStats (s, &stats);
    nsegs = stats.n;
    segs = gv_calloc(2 * nsegs, sizeof(int));

    state.n = 0;
    state.edges = segs;
    gts_surface_foreach_edge(s, addEdge, &state);

    gts_surface_foreach_face(s, cntFace, &nfaces);

    int *faces = gv_calloc(3 * nfaces, sizeof(int));
    int *neigh = gv_calloc(3 * nfaces, sizeof(int));

    statf.faces = faces;
    statf.neigh = neigh;
    gts_surface_foreach_face(s, addFace, &statf);

    sf->nedges = nsegs;
    sf->edges = segs;
    sf->nfaces = nfaces;
    sf->faces = faces;
    sf->neigh = neigh;

    gts_object_destroy (GTS_OBJECT (s));

    return sf;
}

/* Given n points whose coordinates are stored as (x[2*i],x[2*i+1]),
 * compute a Delaunay triangulation of the points.
 * The number of triangles in the triangulation is returned in tris.
 * The return value t is an array of 3*(*tris) integers,
 * with triangle i having points whose indices are t[3*i], t[3*i+1] and t[3*i+2].
 */
int* 
get_triangles (double *x, int n, int* tris)
{
    GtsSurface* s;
    int nfaces = 0;
    fstate statf;

    if (n <= 2) return NULL;

    s = tri(x, NULL, n, NULL, 0, 0);
    if (!s) return NULL;

    gts_surface_foreach_face(s, cntFace, &nfaces);
    statf.faces = gv_calloc(3 * nfaces, sizeof(int));
    gts_surface_foreach_face(s, addTri, &statf);

    gts_object_destroy (GTS_OBJECT (s));

    *tris = nfaces;
    return statf.faces;
}

void 
freeSurface (surface_t* s)
{
    free (s->edges);
    free (s->faces);
    free (s->neigh);
    free(s);
}
#else
static char* err = "Graphviz built without any triangulation library\n";
int* get_triangles (double *x, int n, int* tris)
{
    (void)x;
    (void)n;
    (void)tris;

    agerrorf("get_triangles: %s\n", err);
    return 0;
}
static v_data *delaunay_triangulation(double *x, double *y, int n) {
    (void)x;
    (void)y;
    (void)n;

    agerrorf("delaunay_triangulation: %s\n", err);
    return 0;
}
int *delaunay_tri(double *x, double *y, int n, int* nedges)
{
    (void)x;
    (void)y;
    (void)n;
    (void)nedges;

    agerrorf("delaunay_tri: %s\n", err);
    return 0;
}
surface_t* 
mkSurface (double *x, double *y, int n, int* segs, int nsegs)
{
    (void)x;
    (void)y;
    (void)n;
    (void)segs;
    (void)nsegs;

    agerrorf("mkSurface: %s\n", err);
    return 0;
}
void 
freeSurface (surface_t* s)
{
    (void)s;

    agerrorf("freeSurface: %s\n", err);
}
#endif

static void remove_edge(v_data * graph, int source, int dest)
{
    int i;
    for (i = 1; i < graph[source].nedges; i++) {
	if (graph[source].edges[i] == dest) {
	    graph[source].edges[i] = graph[source].edges[--graph[source].nedges];
	    break;
	}
    }
}

v_data *UG_graph(double *x, double *y, int n) {
    v_data *delaunay;
    int i;
    double dist_ij, dist_ik, dist_jk, x_i, y_i, x_j, y_j;
    int j, k, neighbor_j, neighbor_k;

    if (n == 2) {
	int *edges = gv_calloc(4, sizeof(int));
	delaunay = gv_calloc(n, sizeof(v_data));
	delaunay[0].ewgts = NULL;
	delaunay[0].edges = edges;
	delaunay[0].nedges = 2;
	delaunay[0].edges[0] = 0;
	delaunay[0].edges[1] = 1;
	delaunay[1].edges = edges + 2;
	delaunay[1].ewgts = NULL;
	delaunay[1].nedges = 2;
	delaunay[1].edges[0] = 1;
	delaunay[1].edges[1] = 0;
	return delaunay;
    } else if (n == 1) {
	int *edges = gv_calloc(1, sizeof(int));
	delaunay = gv_calloc(n, sizeof(v_data));
	delaunay[0].ewgts = NULL;
	delaunay[0].edges = edges;
	delaunay[0].nedges = 1;
	delaunay[0].edges[0] = 0;
	return delaunay;
    }

    delaunay = delaunay_triangulation(x, y, n);

    // remove all edges v-u if there is w, neighbor of u or v, that is closer to both u and v than dist(u,v)
    for (i = 0; i < n; i++) {
        x_i = x[i];
        y_i = y[i];
        for (j = 1; j < delaunay[i].nedges;) {
            neighbor_j = delaunay[i].edges[j];
            x_j = x[neighbor_j];
            y_j = y[neighbor_j];
            dist_ij = (x_j - x_i) * (x_j - x_i) + (y_j - y_i) * (y_j - y_i);
            // now look at i'th neighbors to see whether there is a node in the "forbidden region"
            // we will also go through neighbor_j's neighbors when we traverse the edge from its other side
            bool removed = false;
            for (k = 1; k < delaunay[i].nedges && !removed; k++) {
                neighbor_k = delaunay[i].edges[k];
                dist_ik = (x[neighbor_k] - x_i) * (x[neighbor_k] - x_i) +
                    (y[neighbor_k] - y_i) * (y[neighbor_k] - y_i);
                if (dist_ik < dist_ij) {
                    dist_jk = (x[neighbor_k] - x_j) * (x[neighbor_k] - x_j) +
                        (y[neighbor_k] - y_j) * (y[neighbor_k] - y_j);
                    if (dist_jk < dist_ij) {
                        // remove the edge between i and neighbor j
                        delaunay[i].edges[j] = delaunay[i].edges[--delaunay[i].nedges];
                        remove_edge(delaunay, neighbor_j, i);
                        removed = true;
                    }
                }
            }
            if (!removed) {
                j++;
            }
        }
    }
    return delaunay;
}

void freeGraph (v_data * graph)
{
    if (graph != NULL) {
	free(graph[0].edges);
	free(graph[0].ewgts);
	free(graph);
    }
}

void freeGraphData(vtx_data * graph)
{
    if (graph != NULL) {
	free(graph[0].edges);
	free(graph[0].ewgts);
#ifdef DIGCOLA
	free(graph[0].edists);
#endif
	free(graph);
    }
}