summaryrefslogtreecommitdiff
path: root/src/fortune
diff options
context:
space:
mode:
Diffstat (limited to 'src/fortune')
-rw-r--r--src/fortune/defs.h134
-rw-r--r--src/fortune/edgelist.c136
-rw-r--r--src/fortune/geometry.c184
-rw-r--r--src/fortune/heap.c94
-rw-r--r--src/fortune/main.c284
-rw-r--r--src/fortune/memory.c44
-rw-r--r--src/fortune/output.c46
-rw-r--r--src/fortune/voronoi.c103
8 files changed, 1025 insertions, 0 deletions
diff --git a/src/fortune/defs.h b/src/fortune/defs.h
new file mode 100644
index 0000000..b186581
--- /dev/null
+++ b/src/fortune/defs.h
@@ -0,0 +1,134 @@
+#ifndef NULL
+#define NULL 0
+#endif
+#define DELETED -2
+#include <stdlib.h>
+#include <limits.h>
+#include <stdbool.h>
+
+extern int triangulate, sorted, plot, debug;
+
+struct Freenode {
+ struct Freenode *nextfree;
+};
+struct Freelist {
+ struct Freenode *head;
+ int nodesize;
+};
+char *getfree();
+char *myalloc();
+
+extern double xmin, xmax, ymin, ymax, deltax, deltay;
+
+struct Point {
+ double x, y;
+};
+
+/* structure used both for sites and for vertices */
+struct Site {
+ struct Point coord;
+ int sitenbr;
+ int refcnt;
+};
+
+extern struct Site *sites;
+extern int nsites;
+extern int siteidx;
+extern int sqrt_nsites;
+extern int nvertices;
+extern struct Freelist sfl;
+extern struct Site *bottomsite;
+
+char **alloclist;
+int allocnum;
+
+struct Edge {
+ double a, b, c;
+ struct Site *ep[2];
+ struct Site *reg[2];
+ int edgenbr;
+};
+#define le 0
+#define re 1
+extern int nedges;
+extern struct Freelist efl;
+
+int has_endpoint(), right_of();
+struct Site *intersect();
+double dist();
+struct Point PQ_min();
+struct Halfedge *PQextractmin();
+struct Edge *bisect();
+
+struct Halfedge {
+ struct Halfedge *ELleft, *ELright;
+ struct Edge *ELedge;
+ int ELrefcnt;
+ char ELpm;
+ struct Site *vertex;
+ double ystar;
+ struct Halfedge *PQnext;
+};
+
+extern struct Freelist hfl;
+extern struct Halfedge *ELleftend, *ELrightend;
+extern int ELhashsize;
+extern struct Halfedge **ELhash;
+struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
+struct Site *leftreg(), *rightreg();
+
+extern int PQhashsize;
+extern struct Halfedge *PQhash;
+extern struct Halfedge *PQfind();
+extern int PQcount;
+extern int PQmin;
+int PQempty();
+
+extern double *site_list;
+extern double *vert_list;
+extern double *line_list;
+extern unsigned int *edge_list;
+extern unsigned int *dual_list;
+
+extern int vert_count;
+extern int edge_count;
+extern int line_count;
+extern int dual_count;
+extern int site_count;
+
+void makefree(struct Freenode *curr, struct Freelist *fl);
+void voronoi(int triangulate, struct Site *(*nextsite)());
+char *getfree(struct Freelist *fl);
+char *myalloc(unsigned n);
+void freeinit(struct Freelist *fl, int size);
+void ELinitialize();
+struct Halfedge *HEcreate(struct Edge *e, int pm);
+void ELinsert(struct Halfedge *lb, struct Halfedge *new);
+struct Halfedge *ELgethash(int b);
+struct Halfedge *ELleftbnd(struct Point *p);
+void ELdelete(struct Halfedge *he);
+struct Halfedge *ELright(struct Halfedge *he);
+struct Halfedge *ELleft(struct Halfedge *he);
+struct Site *leftreg(struct Halfedge *he);
+struct Site *rightreg(struct Halfedge *he);
+void geominit();
+struct Edge *bisect(struct Site *s1, struct Site *s2);
+struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2);
+int right_of(struct Halfedge *el, struct Point *p);
+void endpoint(struct Edge *e, int lr, struct Site *s);
+double dist(struct Site *s, struct Site *t);
+void makevertex(struct Site *v);
+void deref(struct Site *v);
+void ref(struct Site *v);
+void PQinsert(struct Halfedge *he, struct Site *v, double offset);
+void PQdelete(struct Halfedge *he);
+int PQbucket(struct Halfedge *he);
+int PQempty();
+struct Point PQ_min();
+struct Halfedge *PQextractmin();
+void PQinitialize();
+void out_bisector(struct Edge *e);
+void out_ep(struct Edge *e);
+void out_vertex(struct Site *v);
+void out_site(struct Site *s);
+void out_triple(struct Site *s1, struct Site *s2, struct Site *s3);
diff --git a/src/fortune/edgelist.c b/src/fortune/edgelist.c
new file mode 100644
index 0000000..dfd85a7
--- /dev/null
+++ b/src/fortune/edgelist.c
@@ -0,0 +1,136 @@
+#
+#include "defs.h"
+int nedges;
+struct Freelist efl;
+struct Freelist hfl;
+struct Halfedge *ELleftend, *ELrightend;
+int ELhashsize;
+struct Halfedge **ELhash;
+
+int ntry, totalsearch;
+
+void ELinitialize() {
+ int i;
+ freeinit(&hfl, sizeof **ELhash);
+ ELhashsize = 2 * sqrt_nsites;
+ ELhash = (struct Halfedge **)myalloc(sizeof *ELhash * ELhashsize);
+ for (i = 0; i < ELhashsize; i += 1)
+ ELhash[i] = (struct Halfedge *)NULL;
+ ELleftend = HEcreate((struct Edge *)NULL, 0);
+ ELrightend = HEcreate((struct Edge *)NULL, 0);
+ ELleftend->ELleft = (struct Halfedge *)NULL;
+ ELleftend->ELright = ELrightend;
+ ELrightend->ELleft = ELleftend;
+ ELrightend->ELright = (struct Halfedge *)NULL;
+ ELhash[0] = ELleftend;
+ ELhash[ELhashsize - 1] = ELrightend;
+}
+
+struct Halfedge *HEcreate(e, pm) struct Edge *e;
+int pm;
+{
+ struct Halfedge *answer;
+ answer = (struct Halfedge *)getfree(&hfl);
+ answer->ELedge = e;
+ answer->ELpm = pm;
+ answer->PQnext = (struct Halfedge *)NULL;
+ answer->vertex = (struct Site *)NULL;
+ answer->ELrefcnt = 0;
+ return (answer);
+}
+
+void ELinsert(struct Halfedge *lb, struct Halfedge *new) {
+ new->ELleft = lb;
+ new->ELright = lb->ELright;
+ (lb->ELright)->ELleft = new;
+ lb->ELright = new;
+}
+
+/* Get entry from hash table, pruning any deleted nodes */
+struct Halfedge *ELgethash(b) int b;
+{
+ struct Halfedge *he;
+
+ if (b < 0 || b >= ELhashsize)
+ return ((struct Halfedge *)NULL);
+ he = ELhash[b];
+ if (he == (struct Halfedge *)NULL || he->ELedge != (struct Edge *)DELETED)
+ return (he);
+
+ /* Hash table points to deleted half edge. Patch as necessary. */
+ ELhash[b] = (struct Halfedge *)NULL;
+ if ((he->ELrefcnt -= 1) == 0)
+ makefree((struct Freenode *)he, &hfl);
+ return ((struct Halfedge *)NULL);
+}
+
+struct Halfedge *ELleftbnd(p) struct Point *p;
+{
+ int i, bucket;
+ struct Halfedge *he;
+
+ /* Use hash table to get close to desired halfedge */
+ bucket = (p->x - xmin) / deltax * ELhashsize;
+ if (bucket < 0)
+ bucket = 0;
+ if (bucket >= ELhashsize)
+ bucket = ELhashsize - 1;
+ he = ELgethash(bucket);
+ if (he == (struct Halfedge *)NULL) {
+ for (i = 1; 1; i += 1) {
+ if ((he = ELgethash(bucket - i)) != (struct Halfedge *)NULL)
+ break;
+ if ((he = ELgethash(bucket + i)) != (struct Halfedge *)NULL)
+ break;
+ };
+ totalsearch += i;
+ };
+ ntry += 1;
+ /* Now search linear list of halfedges for the corect one */
+ if (he == ELleftend || (he != ELrightend && right_of(he, p))) {
+ do {
+ he = he->ELright;
+ } while (he != ELrightend && right_of(he, p));
+ he = he->ELleft;
+ } else
+ do {
+ he = he->ELleft;
+ } while (he != ELleftend && !right_of(he, p));
+
+ /* Update hash table and reference counts */
+ if (bucket > 0 && bucket < ELhashsize - 1) {
+ if (ELhash[bucket] != (struct Halfedge *)NULL)
+ ELhash[bucket]->ELrefcnt -= 1;
+ ELhash[bucket] = he;
+ ELhash[bucket]->ELrefcnt += 1;
+ };
+ return (he);
+}
+
+/* This delete routine can't reclaim node, since pointers from hash
+ table may be present. */
+void ELdelete(struct Halfedge *he) {
+ (he->ELleft)->ELright = he->ELright;
+ (he->ELright)->ELleft = he->ELleft;
+ he->ELedge = (struct Edge *)DELETED;
+}
+
+struct Halfedge *ELright(he) struct Halfedge *he;
+{ return (he->ELright); }
+
+struct Halfedge *ELleft(he) struct Halfedge *he;
+{ return (he->ELleft); }
+
+struct Site *leftreg(he) struct Halfedge *he;
+{
+ if (he->ELedge == (struct Edge *)NULL)
+ return (bottomsite);
+ return (he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]);
+}
+
+struct Site *rightreg(he) struct Halfedge *he;
+{
+ if (he->ELedge == (struct Edge *)NULL)
+ return (bottomsite);
+ return (he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]);
+}
diff --git a/src/fortune/geometry.c b/src/fortune/geometry.c
new file mode 100644
index 0000000..131cacf
--- /dev/null
+++ b/src/fortune/geometry.c
@@ -0,0 +1,184 @@
+#
+#include "defs.h"
+#include <math.h>
+
+void geominit() {
+ struct Edge e;
+ double sn;
+
+ freeinit(&efl, sizeof e);
+ nvertices = 0;
+ nedges = 0;
+ sn = nsites + 4;
+ sqrt_nsites = sqrt(sn);
+ deltay = ymax - ymin;
+ deltax = xmax - xmin;
+}
+
+struct Edge *bisect(s1, s2) struct Site *s1, *s2;
+{
+ double dx, dy, adx, ady;
+ struct Edge *newedge;
+
+ newedge = (struct Edge *)getfree(&efl);
+
+ newedge->reg[0] = s1;
+ newedge->reg[1] = s2;
+ ref(s1);
+ ref(s2);
+ newedge->ep[0] = (struct Site *)NULL;
+ newedge->ep[1] = (struct Site *)NULL;
+
+ dx = s2->coord.x - s1->coord.x;
+ dy = s2->coord.y - s1->coord.y;
+ adx = dx > 0 ? dx : -dx;
+ ady = dy > 0 ? dy : -dy;
+ newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+ if (adx > ady) {
+ newedge->a = 1.0;
+ newedge->b = dy / dx;
+ newedge->c /= dx;
+ } else {
+ newedge->b = 1.0;
+ newedge->a = dx / dy;
+ newedge->c /= dy;
+ };
+
+ newedge->edgenbr = nedges;
+ out_bisector(newedge);
+ nedges += 1;
+ return (newedge);
+}
+
+struct Site *intersect(el1, el2) struct Halfedge *el1, *el2;
+{
+ struct Edge *e1, *e2, *e;
+ struct Halfedge *el;
+ double d, xint, yint;
+ int right_of_site;
+ struct Site *v;
+
+ e1 = el1->ELedge;
+ e2 = el2->ELedge;
+ if (e1 == (struct Edge *)NULL || e2 == (struct Edge *)NULL)
+ return ((struct Site *)NULL);
+ if (e1->reg[1] == e2->reg[1])
+ return ((struct Site *)NULL);
+
+ d = e1->a * e2->b - e1->b * e2->a;
+ /* printf("intersect: d=%g\n", d); */
+ if (-1.0e-20 < d && d < 1.0e-20) {
+ return ((struct Site *)NULL);
+ };
+
+ xint = (e1->c * e2->b - e2->c * e1->b) / d;
+ yint = (e2->c * e1->a - e1->c * e2->a) / d;
+
+ if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
+ (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
+ e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
+ el = el1;
+ e = e1;
+ } else {
+ el = el2;
+ e = e2;
+ };
+ right_of_site = xint >= e->reg[1]->coord.x;
+ if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re))
+ return ((struct Site *)NULL);
+
+ v = (struct Site *)getfree(&sfl);
+ v->refcnt = 0;
+ v->coord.x = xint;
+ v->coord.y = yint;
+ return (v);
+}
+
+/* returns 1 if p is to right of halfedge e */
+int right_of(el, p) struct Halfedge *el;
+struct Point *p;
+{
+ struct Edge *e;
+ struct Site *topsite;
+ int right_of_site, above, fast;
+ double dxp, dyp, dxs, t1, t2, t3, yl;
+
+ e = el->ELedge;
+ topsite = e->reg[1];
+ right_of_site = p->x > topsite->coord.x;
+ if (right_of_site && el->ELpm == le)
+ return (1);
+ if (!right_of_site && el->ELpm == re)
+ return (0);
+
+ if (e->a == 1.0) {
+ dyp = p->y - topsite->coord.y;
+ dxp = p->x - topsite->coord.x;
+ fast = 0;
+ if ((!right_of_site && e->b < 0.0) | (right_of_site && e->b >= 0.0)) {
+ above = dyp >= e->b * dxp;
+ fast = above;
+ } else {
+ above = p->x + p->y * e->b > e->c;
+ if (e->b < 0.0)
+ above = !above;
+ if (!above)
+ fast = 1;
+ };
+ if (!fast) {
+ dxs = topsite->coord.x - (e->reg[0])->coord.x;
+ above = e->b * (dxp * dxp - dyp * dyp) <
+ dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
+ if (e->b < 0.0)
+ above = !above;
+ };
+ } else /*e->b==1.0 */
+ {
+ yl = e->c - e->a * p->x;
+ t1 = p->y - yl;
+ t2 = p->x - topsite->coord.x;
+ t3 = yl - topsite->coord.y;
+ above = t1 * t1 > t2 * t2 + t3 * t3;
+ };
+ return (el->ELpm == le ? above : !above);
+}
+
+void endpoint(e, lr, s) struct Edge *e;
+int lr;
+struct Site *s;
+{
+ e->ep[lr] = s;
+ ref(s);
+ if (e->ep[re - lr] == (struct Site *)NULL) {
+ } else {
+ out_ep(e);
+ deref(e->reg[le]);
+ deref(e->reg[re]);
+ makefree((struct Freenode *)e, &efl);
+ }
+}
+
+double dist(s, t) struct Site *s, *t;
+{
+ double dx, dy;
+ dx = s->coord.x - t->coord.x;
+ dy = s->coord.y - t->coord.y;
+ return (sqrt(dx * dx + dy * dy));
+}
+
+void makevertex(v) struct Site *v;
+{
+ v->sitenbr = nvertices;
+ nvertices += 1;
+ out_vertex(v);
+}
+
+void deref(v) struct Site *v;
+{
+ v->refcnt -= 1;
+ if (v->refcnt == 0)
+ makefree((struct Freenode *)v, &sfl);
+}
+
+void ref(v) struct Site *v;
+{ v->refcnt += 1; }
diff --git a/src/fortune/heap.c b/src/fortune/heap.c
new file mode 100644
index 0000000..2ebbd58
--- /dev/null
+++ b/src/fortune/heap.c
@@ -0,0 +1,94 @@
+#
+#include "defs.h"
+int PQhashsize;
+struct Halfedge *PQhash;
+struct Halfedge *PQfind();
+int PQcount;
+int PQmin;
+int PQempty();
+
+void PQinsert(he, v, offset) struct Halfedge *he;
+struct Site *v;
+double offset;
+{
+ struct Halfedge *last, *next;
+ he->vertex = v;
+ ref(v);
+ he->ystar = v->coord.y + offset;
+ last = &PQhash[PQbucket(he)];
+ while ((next = last->PQnext) != (struct Halfedge *)NULL &&
+ (he->ystar > next->ystar ||
+ (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
+ last = next;
+ };
+ he->PQnext = last->PQnext;
+ last->PQnext = he;
+ PQcount += 1;
+}
+
+void PQdelete(he) struct Halfedge *he;
+{
+ struct Halfedge *last;
+
+ if (he->vertex != (struct Site *)NULL) {
+ last = &PQhash[PQbucket(he)];
+ while (last->PQnext != he)
+ last = last->PQnext;
+ last->PQnext = he->PQnext;
+ PQcount -= 1;
+ deref(he->vertex);
+ he->vertex = (struct Site *)NULL;
+ };
+}
+
+int PQbucket(he) struct Halfedge *he;
+{
+ int bucket;
+
+ if (he->ystar < ymin)
+ bucket = 0;
+ else if (he->ystar >= ymax)
+ bucket = PQhashsize - 1;
+ else
+ bucket = (he->ystar - ymin) / deltay * PQhashsize;
+ if (bucket < 0)
+ bucket = 0;
+ if (bucket >= PQhashsize)
+ bucket = PQhashsize - 1;
+ if (bucket < PQmin)
+ PQmin = bucket;
+ return (bucket);
+}
+
+int PQempty() { return (PQcount == 0); }
+
+struct Point PQ_min() {
+ struct Point answer;
+
+ while (PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {
+ PQmin += 1;
+ };
+ answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
+ answer.y = PQhash[PQmin].PQnext->ystar;
+ return (answer);
+}
+
+struct Halfedge *PQextractmin() {
+ struct Halfedge *curr;
+ curr = PQhash[PQmin].PQnext;
+ PQhash[PQmin].PQnext = curr->PQnext;
+ PQcount -= 1;
+ return (curr);
+}
+
+void PQinitialize() {
+ int i;
+// struct Point *s;
+
+ PQcount = 0;
+ PQmin = 0;
+ PQhashsize = 4 * sqrt_nsites;
+ PQhash = (struct Halfedge *)myalloc(PQhashsize * sizeof *PQhash);
+ for (i = 0; i < PQhashsize; i += 1)
+ PQhash[i].PQnext = (struct Halfedge *)NULL;
+}
diff --git a/src/fortune/main.c b/src/fortune/main.c
new file mode 100644
index 0000000..7571945
--- /dev/null
+++ b/src/fortune/main.c
@@ -0,0 +1,284 @@
+#
+#include <stdio.h>
+#include <stdint.h>
+#include "defs.h"
+#include <math.h>
+struct Site *nextone();
+
+int triangulate, sorted, plot, debug;
+double xmin, xmax, ymin, ymax, deltax, deltay;
+int vert_count, edge_count, line_count, dual_count, site_count;
+double *site_list, *vert_list, *line_list;
+unsigned int *edge_list, *dual_list;
+
+/* sort sites on y, then x, coord */
+int scomp(s1, s2) struct Point *s1, *s2;
+{
+ if (s1->y < s2->y)
+ return (-1);
+ if (s1->y > s2->y)
+ return (1);
+ if (s1->x < s2->x)
+ return (-1);
+ if (s1->x > s2->x)
+ return (1);
+ return (0);
+}
+
+/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
+void readsites(double *lattice) {
+
+ sites = (struct Site *)myalloc(nsites * sizeof *sites);
+ for (unsigned int j = 0; j < nsites; j++) {
+ sites[j].coord.x = lattice[2 * j];
+ sites[j].coord.y = lattice[2 * j + 1];
+ sites[j].sitenbr = j;
+ sites[j].refcnt = 0;
+ }
+ qsort(sites, nsites, sizeof(struct Site), scomp);
+ xmin = sites[0].coord.x;
+ xmax = sites[0].coord.x;
+ for (unsigned int i = 1; i < nsites; i += 1) {
+ if (sites[i].coord.x < xmin)
+ xmin = sites[i].coord.x;
+ if (sites[i].coord.x > xmax)
+ xmax = sites[i].coord.x;
+ };
+ ymin = sites[0].coord.y;
+ ymax = sites[nsites - 1].coord.y;
+}
+
+unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) {
+ unsigned int ndup = 0;
+ bool *duplicates = (bool *)calloc(ne, sizeof(bool));
+ for (unsigned int i = 0; i < ne; i++) {
+ unsigned int v1 = etv[2 * i];
+ unsigned int v2 = etv[2 * i + 1];
+ for (unsigned int j = 0; j < i; j++) {
+ unsigned int w1 = etv[2 * j];
+ unsigned int w2 = etv[2 * j + 1];
+ if (w1 == v1 && w2 == v2) {
+ duplicates[j] = true;
+ ndup++;
+ break;
+ }
+ }
+ }
+ unsigned int newsize = (int)ne - (int)ndup;
+ unsigned int count = 0;
+ for (unsigned int i = 0; i < ne; i++) {
+ if (!duplicates[i]) {
+ etv[2 * count] = etv[2 * i];
+ etv[2 * count + 1] = etv[2 * i + 1];
+ count++;
+ }
+ }
+ free(duplicates);
+ return newsize;
+}
+
+intptr_t *run_voronoi(unsigned int num, double *lattice, double xmin, double xmax, double ymin, double ymax)
+{
+ struct Site *(*next)();
+
+ sorted = 0;
+ triangulate = 0;
+ plot = 0;
+ debug = 0;
+
+ alloclist = (char **)malloc(9 * num * sizeof(char *));
+ allocnum = 0;
+
+ freeinit(&sfl, sizeof *sites);
+
+ unsigned int eff_num = 9 * num;
+ double *eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
+
+ for (unsigned int i = 0; i < num; i++) {
+ // original sites - our baby boys
+ eff_lattice[2*i] = lattice[2*i];
+ eff_lattice[2*i+1] = lattice[2*i+1];
+
+ // sites to the right
+ eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax;
+ eff_lattice[2*(num+i)] = lattice[2*i];
+
+ // sites to the left
+ eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+ eff_lattice[2*(2*num+i)] = lattice[2*i];
+
+ // sites to the top
+ eff_lattice[2*(3*num+i)+1] = lattice[2*i+1];
+ eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax;
+
+ // sites to the bottom
+ eff_lattice[2*(4*num+i)+1] = lattice[2*i+1];
+ eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin;
+
+ // sites to the upper right
+ eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax;
+ eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax;
+
+ // sites to the upper left
+ eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+ eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax;
+
+ // sites to the lower left
+ eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+ eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax;
+
+ // sites to the lower right
+ eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin;
+ eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax;
+ }
+
+ nsites = eff_num;
+ readsites(eff_lattice);
+ free(eff_lattice);
+ next = nextone;
+
+ siteidx = 0;
+ geominit();
+
+ site_list = malloc(2 * nsites * sizeof(double));
+ vert_list = NULL;
+ line_list = NULL;
+ edge_list = NULL;
+ dual_list = NULL;
+ site_count = 0;
+ vert_count = 0;
+ edge_count = 0;
+ line_count = 0;
+ dual_count = 0;
+
+ voronoi(triangulate, next);
+
+ unsigned int real_vert_count = vert_count;
+ unsigned int real_edge_count = edge_count;
+ unsigned int real_dual_count = dual_count;
+ double *real_vert_list = vert_list;
+ unsigned int *real_edge_list = edge_list;
+ unsigned int *real_dual_list = dual_list;
+
+ {
+ real_vert_count = 0;
+ real_edge_count = 0;
+ real_dual_count = 0;
+ real_vert_list = (double *)malloc(2 * vert_count * sizeof(double));
+ real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int));
+ real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int));
+ unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int));
+ unsigned int q = 0;
+ int *new_vert_ind = (int *)malloc(vert_count * sizeof(int));
+
+ for (unsigned int i = 0; i < vert_count; i++) {
+ unsigned int t1, t2, t3;
+ t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2];
+ if (t1 >= num && t2 >= num && t3 >= num) {
+ new_vert_ind[i] = -1;
+ }
+ else if (t1 < num && t2 < num && t3 < num) {
+ real_vert_list[2*real_vert_count] = vert_list[2*i];
+ real_vert_list[2*real_vert_count+1] = vert_list[2*i+1];
+ real_dual_list[3*real_vert_count] = t1;
+ real_dual_list[3*real_vert_count+1] = t2;
+ real_dual_list[3*real_vert_count+2] = t3;
+ new_vert_ind[i] = real_vert_count;
+ real_vert_count++;
+ }
+ else {
+ unsigned int tt1, tt2, tt3;
+ tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num;
+ unsigned int s1, s2, s3;
+ s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3);
+ s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3);
+ s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2));
+
+ bool matched = false;
+ unsigned int ass_vert;
+
+ for (unsigned int j = 0; j < q; j++) {
+ if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) {
+ matched = true;
+ ass_vert = triangle_analyzed[4*j+3];
+ break;
+ }
+ }
+
+ if (matched) {
+ new_vert_ind[i] = ass_vert;
+ } else {
+ triangle_analyzed[4*q] = s1;
+ triangle_analyzed[4*q+1] = s2;
+ triangle_analyzed[4*q+2] = s3;
+ triangle_analyzed[4*q+3] = real_vert_count;
+ q++;
+ real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax);
+ real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax);
+ real_dual_list[3*real_vert_count] = t1 % num;
+ real_dual_list[3*real_vert_count+1] = t2 % num;
+ real_dual_list[3*real_vert_count+2] = t3 % num;
+ new_vert_ind[i] = real_vert_count;
+ real_vert_count++;
+ }
+ }
+ }
+ for (unsigned int i = 0; i < edge_count; i++) {
+ unsigned int v1, v2;
+ v1 = edge_list[2 * i];
+ v2 = edge_list[2 * i+1];
+ if (v1 != UINT_MAX && v2 != UINT_MAX) {
+ if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) {
+ unsigned int nv1 = new_vert_ind[v1];
+ unsigned int nv2 = new_vert_ind[v2];
+ real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2;
+ real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1;
+ real_edge_count++;
+ }
+ }
+ }
+ unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list);
+ real_edge_count = new_edge_count;
+ free(triangle_analyzed);
+ free(new_vert_ind);
+ free(dual_list);
+ free(edge_list);
+ free(vert_list);
+ }
+
+
+ intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t));
+ output[0] = (intptr_t)malloc(4 * sizeof(unsigned int));
+ ((unsigned int *)output[0])[0] = real_vert_count;
+ ((unsigned int *)output[0])[1] = real_edge_count;
+ ((unsigned int *)output[0])[2] = line_count;
+ ((unsigned int *)output[0])[3] = real_vert_count;
+ output[1] = (intptr_t)site_list;
+ output[2] = (intptr_t)real_vert_list;
+ output[3] = (intptr_t)real_edge_list;
+ output[4] = (intptr_t)line_list;
+ output[5] = (intptr_t)real_dual_list;
+
+ free(ELhash);
+ free(PQhash);
+ free(sites);
+ for (unsigned int i = 0; i < allocnum; i++) {
+ free(alloclist[i]);
+ }
+ free(alloclist);
+
+ return output;
+}
+
+/* return a single in-storage site */
+struct Site *nextone() {
+ for (; siteidx < nsites; siteidx += 1) {
+ if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x ||
+ sites[siteidx].coord.y != sites[siteidx - 1].coord.y) {
+ siteidx += 1;
+ return (&sites[siteidx - 1]);
+ };
+ };
+ return ((struct Site *)NULL);
+}
+
diff --git a/src/fortune/memory.c b/src/fortune/memory.c
new file mode 100644
index 0000000..807d1b0
--- /dev/null
+++ b/src/fortune/memory.c
@@ -0,0 +1,44 @@
+#
+#include "defs.h"
+#include <stdio.h>
+
+void freeinit(struct Freelist *fl, int size) {
+ fl->head = (struct Freenode *)NULL;
+ fl->nodesize = size;
+}
+
+char *getfree(fl) struct Freelist *fl;
+{
+ int i;
+ struct Freenode *t;
+ if (fl->head == (struct Freenode *)NULL) {
+ t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize);
+ alloclist[allocnum] = t;
+ allocnum++;
+ for (i = 0; i < sqrt_nsites; i += 1)
+ makefree((struct Freenode *)((char *)t + i * fl->nodesize), fl);
+ };
+ t = fl->head;
+ fl->head = (fl->head)->nextfree;
+ return ((char *)t);
+}
+
+void makefree(curr, fl) struct Freenode *curr;struct Freelist *fl;
+{
+ curr->nextfree = fl->head;
+ fl->head = curr;
+}
+
+int total_alloc;
+char *myalloc(n) unsigned n;
+{
+ char *t;
+ if ((t = malloc(n)) == (char *)0) {
+ fprintf(stderr,
+ "Insufficient memory processing site %d (%d bytes in use)\n",
+ siteidx, total_alloc);
+ exit(1);
+ };
+ total_alloc += n;
+ return (t);
+}
diff --git a/src/fortune/output.c b/src/fortune/output.c
new file mode 100644
index 0000000..d496feb
--- /dev/null
+++ b/src/fortune/output.c
@@ -0,0 +1,46 @@
+#
+#include "defs.h"
+#include <stdio.h>
+double pxmin, pxmax, pymin, pymax, cradius;
+
+void out_bisector(e) struct Edge *e;
+{
+ if (line_count % nsites == 0) line_list = realloc(line_list, 3 * (nsites + line_count) * sizeof(double));
+ line_list[3*line_count] = e->a;
+ line_list[3*line_count+1] = e->b;
+ line_list[3*line_count+2] = e->c;
+ line_count++;
+}
+
+void out_ep(e) struct Edge *e;
+{
+ if (edge_count % nsites == 0) edge_list = realloc(edge_list, 2 * (nsites + edge_count) * sizeof(unsigned int));
+ edge_list[2 * edge_count] = (e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : UINT_MAX);
+ edge_list[2 * edge_count + 1] = (e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : UINT_MAX);
+ edge_count++;
+}
+
+void out_vertex(v) struct Site *v;
+{
+ if (vert_count % nsites == 0) vert_list = realloc(vert_list, 2 * (nsites + vert_count) * sizeof(double));
+ vert_list[2 * vert_count] = v->coord.x;
+ vert_list[2 * vert_count + 1] = v->coord.y;
+ vert_count++;
+}
+
+void out_site(s) struct Site *s;
+{
+ site_list[2 * site_count] = s->coord.x;
+ site_list[2 * site_count + 1] = s->coord.y;
+ site_count++;
+}
+
+void out_triple(s1, s2, s3) struct Site *s1, *s2, *s3;
+{
+ if (dual_count % nsites == 0) dual_list = realloc(dual_list, 3 * (nsites + dual_count) * sizeof(unsigned int));
+ dual_list[dual_count * 3] = s1->sitenbr;
+ dual_list[dual_count * 3 + 1] = s2->sitenbr;
+ dual_list[dual_count * 3 + 2] = s3->sitenbr;
+ dual_count++;
+}
+
diff --git a/src/fortune/voronoi.c b/src/fortune/voronoi.c
new file mode 100644
index 0000000..dc9945f
--- /dev/null
+++ b/src/fortune/voronoi.c
@@ -0,0 +1,103 @@
+#
+#include "defs.h"
+
+struct Site *sites;
+int nsites;
+int siteidx;
+int sqrt_nsites;
+int nvertices;
+struct Freelist sfl;
+struct Site *bottomsite;
+
+/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
+ deltax, deltay (can all be estimates).
+ Performance suffers if they are wrong; better to make nsites,
+ deltax, and deltay too big than too small. (?) */
+
+void voronoi(triangulate, nextsite) int triangulate;
+struct Site *(*nextsite)();
+{
+ struct Site *newsite, *bot, *top, *temp, *p;
+ struct Site *v;
+ struct Point newintstar;
+ int pm;
+ struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
+ struct Edge *e;
+
+ PQinitialize();
+ bottomsite = (*nextsite)();
+ out_site(bottomsite);
+ ELinitialize();
+
+ newsite = (*nextsite)();
+ while (1) {
+ if (!PQempty())
+ newintstar = PQ_min();
+
+ if (newsite != (struct Site *)NULL &&
+ (PQempty() || newsite->coord.y < newintstar.y ||
+ (newsite->coord.y == newintstar.y &&
+ newsite->coord.x < newintstar.x))) {/* new site is smallest */
+ out_site(newsite);
+ lbnd = ELleftbnd(&(newsite->coord));
+ rbnd = ELright(lbnd);
+ bot = rightreg(lbnd);
+ e = bisect(bot, newsite);
+ bisector = HEcreate(e, le);
+ ELinsert(lbnd, bisector);
+ if ((p = intersect(lbnd, bisector)) != (struct Site *)NULL) {
+ PQdelete(lbnd);
+ PQinsert(lbnd, p, dist(p, newsite));
+ };
+ lbnd = bisector;
+ bisector = HEcreate(e, re);
+ ELinsert(lbnd, bisector);
+ if ((p = intersect(bisector, rbnd)) != (struct Site *)NULL) {
+ PQinsert(bisector, p, dist(p, newsite));
+ };
+ newsite = (*nextsite)();
+ } else if (!PQempty())
+ /* intersection is smallest */
+ {
+ lbnd = PQextractmin();
+ llbnd = ELleft(lbnd);
+ rbnd = ELright(lbnd);
+ rrbnd = ELright(rbnd);
+ bot = leftreg(lbnd);
+ top = rightreg(rbnd);
+ out_triple(bot, top, rightreg(lbnd));
+ v = lbnd->vertex;
+ makevertex(v);
+ endpoint(lbnd->ELedge, lbnd->ELpm, v);
+ endpoint(rbnd->ELedge, rbnd->ELpm, v);
+ ELdelete(lbnd);
+ PQdelete(rbnd);
+ ELdelete(rbnd);
+ pm = le;
+ if (bot->coord.y > top->coord.y) {
+ temp = bot;
+ bot = top;
+ top = temp;
+ pm = re;
+ }
+ e = bisect(bot, top);
+ bisector = HEcreate(e, pm);
+ ELinsert(llbnd, bisector);
+ endpoint(e, re - pm, v);
+ deref(v);
+ if ((p = intersect(llbnd, bisector)) != (struct Site *)NULL) {
+ PQdelete(llbnd);
+ PQinsert(llbnd, p, dist(p, bot));
+ };
+ if ((p = intersect(bisector, rrbnd)) != (struct Site *)NULL) {
+ PQinsert(bisector, p, dist(p, bot));
+ };
+ } else
+ break;
+ };
+
+ for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
+ e = lbnd->ELedge;
+ out_ep(e);
+ };
+}