diff options
Diffstat (limited to 'src/fortune')
| -rw-r--r-- | src/fortune/defs.h | 134 | ||||
| -rw-r--r-- | src/fortune/edgelist.c | 136 | ||||
| -rw-r--r-- | src/fortune/geometry.c | 184 | ||||
| -rw-r--r-- | src/fortune/heap.c | 94 | ||||
| -rw-r--r-- | src/fortune/main.c | 284 | ||||
| -rw-r--r-- | src/fortune/memory.c | 44 | ||||
| -rw-r--r-- | src/fortune/output.c | 46 | ||||
| -rw-r--r-- | src/fortune/voronoi.c | 103 |
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); + }; +} |