diff options
Diffstat (limited to 'test/ccured_olden/voronoi/newvor.c')
| -rw-r--r-- | test/ccured_olden/voronoi/newvor.c | 711 |
1 files changed, 711 insertions, 0 deletions
diff --git a/test/ccured_olden/voronoi/newvor.c b/test/ccured_olden/voronoi/newvor.c new file mode 100644 index 00000000..9b5a4408 --- /dev/null +++ b/test/ccured_olden/voronoi/newvor.c @@ -0,0 +1,711 @@ +/* For copyright information, see olden_v1.0/COPYRIGHT */ + +#ifdef SS_PLAIN +#include "ssplain.h" +#else SS_PLAIN +#include <cm/cmmd.h> +#include <fcntl.h> +#endif SS_PLAIN + +#define DEFINE_GLOBALS +#include "defines.h" +/* WARNING! Don't use LOCAL on QUAD_EDGE as cache addresses are not aligned */ + +int flag; + +QUAD_EDGE connect_left(a, b) +QUAD_EDGE a,b; +{ + VERTEX_PTR t1,t2; + QUAD_EDGE ans,lnexta; + + t1=dest(a); + lnexta=lnext(a); + t2=orig(b); + ans = makeedge(t1,t2/*dest(a), orig(b)*/); + splice(ans, lnexta); + splice(sym(ans), b); + return(ans); +} + +QUAD_EDGE connect_right(a, b) + QUAD_EDGE a,b; +{ + VERTEX_PTR t1,t2; + QUAD_EDGE ans, oprevb; + + t1=dest(a); + t2=orig(b); + + oprevb=oprev(b); + ans = makeedge(t1,t2); + splice(ans, sym(a)); + splice(sym(ans), oprevb); + return(ans); +} + +/****************************************************************/ +/* Top-level Delaunay Triangulation Procedure */ +/****************************************************************/ + +QUAD_EDGE build_delaunay_triangulation(tree,extra) + /* builds delaunay triangulation. + va is an array of vertices, from 0 to size. Each vertex consists of + a vector and a data pointer. edge is a pointer to an + edge on the convex hull of the constructed delaunay triangulation. */ + + VERTEX_PTR tree,extra; +{ + EDGE_PAIR retval; + + retval=build_delaunay(tree,extra); + return retval.left; +} + +VERTEX_PTR get_low(tree) + VERTEX_PTR tree; +{ + VERTEX_PTR temp; + while((temp=tree->left)) tree=temp; + return tree; +} + +/****************************************************************/ +/* Recursive Delaunay Triangulation Procedure */ +/* Contains modifications for axis-switching division. */ +/****************************************************************/ + +EDGE_PAIR build_delaunay(tree,extra) + VERTEX_PTR tree,extra; +{ + QUAD_EDGE a,b,c,ldo,rdi,ldi,rdo; + EDGE_PAIR retval; + VERTEX_PTR maxx, minx; + VERTEX_PTR s1, s2, s3; + + EDGE_PAIR delleft, delright; + + if (tree && tree->right && tree->left) + { + /* more than three elements; do recursion */ + minx = get_low(tree); maxx = extra; + + delright = build_delaunay(tree->right,extra); + delleft = build_delaunay(tree->left, tree); + ldo = delleft.left; ldi=delleft.right; + rdi=delright.left; rdo=delright.right; + + retval=do_merge(ldo, ldi, rdi, rdo); + ldo = retval.left; + rdo = retval.right; + while (orig(ldo) != minx) { ldo = rprev(ldo); } + while (orig(rdo) != maxx) { rdo = lprev(rdo); } + retval.left = ldo; + retval.right = rdo; + } + else if (!tree) + { + printf("ERROR: Only 1 point!\n"); + exit(-1); + } + else if (!tree->left) + { + /* two points */ + a = makeedge(tree, extra); + retval.left = a; + retval.right = sym(a); + } + else + { /* tree->left, !tree->right */ /* three points */ + /* 3 cases: triangles of 2 orientations, and 3 points on a line. */ + s1 = tree->left; + s2 = tree; + s3 = extra; + a = makeedge(s1, s2); + b = makeedge(s2, s3); + splice(sym(a), b); + c = connect_left(b, a); + if (ccw(s1, s3, s2)) + { + retval.left = sym(c); + retval.right = c; + } + else + { + retval.left = a; + retval.right = sym(b); + if (!ccw(s1, s2, s3)) deleteedge(c); /* colinear */ + } + } + return retval; +} + +/****************************************************************/ +/* Quad-edge storage allocation */ +/****************************************************************/ +QUAD_EDGE next_edge, avail_edge; +#ifndef NEW +#define NYL -1 +#else +#define NYL NULL +#endif + +void +free_edge(e) + QUAD_EDGE e; +{ + e = (QUAD_EDGE) ((int) e ^ ((int) e & ANDF)); + onext(e) = avail_edge; + avail_edge = e; +} + +void +deleteedge(e) + /*disconnects e from the rest of the structure and destroys it. */ + QUAD_EDGE e; +{ + QUAD_EDGE f; + f=oprev(e); + splice(e, f); + f=oprev(sym(e)); + splice(sym(e),f); + free_edge(e); +} + +void +delete_all_edges() +{ + next_edge= 0; + avail_edge = NYL; +} + +QUAD_EDGE alloc_edge() +{ + QUAD_EDGE ans; + + if (avail_edge == NYL) + { + ans = (QUAD_EDGE) malloc(4*ALLOC_SIZE); +#ifdef OUT + if ((int) ans & ANDF) +#else !OUT + if (!(int) ans & 0xF) +#endif OUT + { + printf("Aborting in alloc_edge, ans = 0x%x\n",(unsigned)ans); + exit(-1); + } + } + else ans = (QUAD_EDGE) avail_edge, avail_edge = onext(avail_edge); + return ans; +} + +/****************************************************************/ +/* Geometric primitives */ +/****************************************************************/ + +BOOLEAN incircle(a,b,c,d) + /* incircle, as in the Guibas-Stolfi paper. */ + VERTEX_PTR a,b,c,d; +{ + double adx, ady, bdx, bdy, cdx, cdy, dx, dy, anorm, bnorm, cnorm, dnorm; + double dret ; + + dx = X(d); dy = Y(d); dnorm = NORM(d); + adx = X(a) - dx; ady = Y(a) - dy; anorm = NORM(a); + bdx = X(b) - dx; bdy = Y(b) - dy; bnorm = NORM(b); + cdx = X(c) - dx; cdy = Y(c) - dy; cnorm = NORM(c); + dret = (anorm - dnorm) * (bdx * cdy - bdy * cdx); + dret += (bnorm - dnorm) * (cdx * ady - cdy * adx); + dret += (cnorm - dnorm) * (adx * bdy - ady * bdx); + return( (0.0 < dret) ? TRUE : FALSE ); +} + +BOOLEAN ccw(a,b,c) + /* TRUE iff A, B, C form a counterclockwise oriented triangle */ + VERTEX_PTR a,b,c; +{ + double dret ; + double xa,ya,xb,yb,xc,yc; + + xa=X(a); ya=Y(a); + xb=X(b); yb=Y(b); + xc=X(c); yc=Y(c); + + dret = (xa-xc)*(yb-yc) - (xb-xc)*(ya-yc); + return( (dret > 0.0)? TRUE : FALSE); +} + +/****************************************************************/ +/* Quad-edge manipulation primitives */ +/****************************************************************/ +QUAD_EDGE makeedge(origin, destination) + VERTEX_PTR origin, destination; +{ + QUAD_EDGE temp, ans; + temp = alloc_edge(); + ans = temp; + + onext(temp) = ans; + orig(temp) = origin; + temp = (QUAD_EDGE) ((int) temp+SIZE); + onext(temp) = (QUAD_EDGE) ((int) ans + 3*SIZE); + temp = (QUAD_EDGE) ((int) temp+SIZE); + onext(temp) = (QUAD_EDGE) ((int) ans + 2*SIZE); + orig(temp) = destination; + temp = (QUAD_EDGE) ((int) temp+SIZE); + onext(temp) = (QUAD_EDGE) ((int) ans + 1*SIZE); + return(ans); +} + +void +splice(a,b) + QUAD_EDGE a, b; +{ + QUAD_EDGE alpha, beta, temp; + QUAD_EDGE t1; + + alpha = rot(onext(a)); + beta = rot(onext(b)); + + t1 = onext(beta); + temp = onext(alpha); + + onext(alpha) = t1; + onext(beta) = temp; + + temp = onext(a); + t1 = onext(b); + onext(b) = temp; + onext(a) = t1; +} + +void +swapedge(e) + QUAD_EDGE e; +{ + QUAD_EDGE a,b,syme,lnexttmp; + VERTEX_PTR a1,b1; + + a = oprev(e); + syme = sym(e); + b = oprev(syme); + splice(e, a); + splice(syme, b); + lnexttmp = lnext(a); + splice(e, lnexttmp); + lnexttmp = lnext(b); + splice(syme, lnexttmp); + a1 = dest(a); + b1 = dest(b); + orig(e) = a1; + dest(e) = b1; +} + +/****************************************************************/ +/* The Merge Procedure. */ +/****************************************************************/ + +int valid(l,basel) + QUAD_EDGE l,basel; +{ + VERTEX_PTR t1,t2,t3; + + t1=orig(basel); + t3=dest(basel); + + t2=dest(l); + return ccw(t1,t2,t3); +} + +void dump_quad(ptr) + QUAD_EDGE ptr; +{ + int i; + QUAD_EDGE j; + VERTEX_PTR v; + + ptr = (QUAD_EDGE) ((int) ptr & ~ANDF); + chatting("Entered DUMP_QUAD: ptr=0x%x\n",ptr); + for (i=0; i<4; i++) + { + j=onext(((QUAD_EDGE) (ptr+i))); + v = orig(j); + chatting("DUMP_QUAD: ptr=0x%x onext=0x%x,v=0x%x\n",ptr+i,j,v); + } +} + + +EDGE_PAIR do_merge(ldo, ldi, rdi, rdo) + QUAD_EDGE ldi, rdi, ldo, rdo; +{ + int rvalid, lvalid; + QUAD_EDGE basel,lcand,rcand,t; + VERTEX_PTR t1,t2,t3; + + while (TRUE) + { + t3=orig(rdi); + + t1=orig(ldi); + t2=dest(ldi); + + while (ccw(t1,t2,t3/*orig(ldi), dest(ldi), orig(rdi)*/)) + { + ldi = lnext(ldi); + + t1=orig(ldi); + t2=dest(ldi); + } + + t2=dest(rdi); + + if (ccw(t2,t3,t1/*dest(rdi), orig(rdi), orig(ldi)*/)) + { + rdi = rprev(rdi); + } + else + { + break; + } + } + + basel = connect_left(sym(rdi), ldi); + + lcand = rprev(basel); + rcand = oprev(basel); + t1 = orig(basel); + t2 = dest(basel); + + if (t1/*orig(basel)*/ == orig(rdo)) rdo = basel; + if (t2/*dest(basel)*/ == orig(ldo)) ldo = sym(basel); + + while (TRUE) + { + VERTEX_PTR v1,v2,v3,v4; + + /*chatting("valid site 1,lcand=0x%x,basel=0x%x\n",lcand,basel);*/ + /*dump_quad(lcand);*/ + t=onext(lcand); + if (valid(t,basel)) + { + v4=orig(basel); + + v1=dest(lcand); + v3=orig(lcand); + + v2=dest(t); + while (incircle(v1,v2,v3,v4 + /*dest(lcand), dest(t), orig(lcand), orig(basel)*/)) + { + deleteedge(lcand); + lcand = t; + + t = onext(lcand); + v1=dest(lcand); + v3=orig(lcand); + + v2=dest(t); + } + } + + /*chatting("valid site 2\n");*/ + t=oprev(rcand); + if (valid(t,basel)) { + v4=dest(basel); + v1=dest(t); + v2=dest(rcand); + v3=orig(rcand); + while (incircle(v1,v2,v3,v4 + /*dest(t), dest(rcand), orig(rcand), dest(basel)*/)) + { + deleteedge(rcand); + rcand = t; + t = oprev(rcand); + v2=dest(rcand); + v3=orig(rcand); + v1=dest(t); + } + } + /*chatting("Valid sites 3,4\n");*/ + lvalid = valid(lcand,basel); + /*chatting("Valid sites 3,4\n");*/ + rvalid = valid(rcand,basel); + if ((! lvalid) && (! rvalid)) + { + EDGE_PAIR retval; + retval.left = ldo; retval.right = rdo; return retval; + } + v1=dest(lcand); + v2=orig(lcand); + v3=orig(rcand); + v4=dest(rcand); + if (!lvalid || + (rvalid && incircle(v1,v2,v3,v4 + /*dest(lcand), orig(lcand), + orig(rcand), dest(rcand)*/))) + { + basel = connect_left(rcand, sym(basel)); + rcand = lnext(sym(basel)); + } + else + { + basel = sym(connect_right(lcand, basel)); + lcand = rprev(basel); + } + } +} + + + +#define DEFINE_GLOBALS +#define CONST_m1 10000 +#define CONST_b 31415821 +#define RANGE 100 + +struct EDGE_STACK *allocate_stack(); +struct get_point get_points(); +int loop = 0, randum = 1, filein = 0 , fileout = 1, statistics = 0; + +void in_order(tree) + VERTEX_PTR tree; +{ + VERTEX_PTR tleft, tright; + double x, y; + + if (!tree) { + return; + } + + x = X(tree); + y = Y(tree); + chatting("X=%f, Y=%f\n",x, y); + tleft = tree->left; + in_order(tleft); + tright = tree->right; + in_order(tright); +} + +int mult(int p, int q) +{ + int p1, p0, q1, q0; + + p1=p/CONST_m1; p0=p%CONST_m1; + q1=q/CONST_m1; q0=q%CONST_m1; + return (((p0*q1+p1*q0) % CONST_m1)*CONST_m1+p0*q0); +} + +int skiprand(int seed, int n) +/* Generate the nth random # */ +{ + for (; n; n--) +#ifdef SS_PLAIN + seed = olden_random(seed); +#else SS_PLAIN + seed=random(seed); +#endif SS_PLAIN + return seed; +} + + +#ifdef SS_PLAIN +int olden_random(int seed) +#else SS_PLAIN +int random(int seed) +#endif SS_PLAIN +{ + seed = (mult(seed,CONST_b)+1); + return seed; +} + + +void +print_extra(extra) + VERTEX_PTR extra; +{ + + double x, y; + x = X(extra); + y = Y(extra); + chatting("X=%f, Y=%f\n",x, y); +} + + +void +main(argc,argv) + int argc; + char *argv[]; +{ + struct EDGE_STACK *my_stack; + struct get_point point, extra; + QUAD_EDGE edge; + int n, retained; + to_lincoln = to_off = to_3d_out = to_color = 0; + voronoi = delaunay = 1; interactive = ahost = 0 ; + retained = 0; + + filestuff(); + chatting("argc = %d\n",argc); + n = dealwithargs(argc, argv); + + chatting("getting %d points\n", n); + extra=get_points(1,1.0,n,1023); + point=get_points(n-1,extra.curmax,n-1,extra.seed); + chatting("Done getting points\n"); + num_vertices = n + 1; + my_stack=allocate_stack(num_vertices); + +//#ifdef SS_PLAIN +// ssplain_alloc_stats(); +//#endif SS_PLAIN + + if (flag) in_order(point.v); + if (flag) print_extra(extra.v); + chatting("Doing voronoi on %d nodes\n", n); + + edge=build_delaunay_triangulation(point.v,extra.v); + + chatting("Elapsed time %f\n", CMMD_node_timer_elapsed(0)); + if (flag) output_voronoi_diagram(edge,n,my_stack); + + exit(0); +} + +struct EDGE_STACK *allocate_stack(num_vertices) + int num_vertices; +{ + struct EDGE_STACK *my_stack; + + num_edgeparts = 12*num_vertices; + my_alloc(my_stack, struct EDGE_STACK, 1); + my_alloc(my_stack->elts, QUAD_EDGE , num_edgeparts); + my_stack->stack_size = num_edgeparts/2; + return my_stack; +} + +void +free_all(cont,my_stack) + int cont; + struct EDGE_STACK *my_stack; +{ + free(my_stack->elts); + free(my_stack); +} + +struct get_point get_points(n,curmax,i,seed) + int n; + double curmax; + int i,seed; +{ + VERTEX_PTR node; + + struct get_point point; + + if (n<1 || i<=n/2) { + point.v = NULL; + point.curmax=curmax; + point.seed = seed; + return point; + } +#ifdef OUT + chatting("Get points: %d, %f, %d\n",n,curmax,i); +#endif OUT + + point = get_points(n/2,curmax,i,seed); + /*chatting("rec call n=%d\n",n);*/ + i -= n/2; + + node = (VERTEX_PTR) mymalloc(sizeof(struct VERTEX)); + + /*chatting("Get points past alloc,n=%d\n",n);*/ + X(node) = point.curmax * exp(log(drand(point.seed))/i); + Y(node) = drand(point.seed); + NORM(node) = X(node)*X(node) + Y(node)*Y(node); + node->right = point.v; + /*chatting("node = 0x%x\n",node);*/ + point = get_points(n/2,X(node),i-1,point.seed); + + node->left = point.v; + point.v = node; + return point; +} + + +/****************************************************************/ +/* Voronoi Diagram Routines. */ +/****************************************************************/ + +/****************************************************************/ +/* Graph Traversal Routines */ +/****************************************************************/ + +QUAD_EDGE pop_edge(x) + struct EDGE_STACK *x; +{ + int a=x->ptr--; + return (x)->elts[a]; +} + +void +push_edge(stack,edge) + struct EDGE_STACK *stack; + QUAD_EDGE edge; +{ + int a; + /*chatting("pushing edge \n");*/ + if (stack->ptr == stack->stack_size) { + printf("cannot push onto stack: stack is too large\n"); + } + else { + a = stack->ptr; + a++; + stack->ptr = a; + stack->elts[a] = edge; + } +} + +void +push_ring(stack, edge) + struct EDGE_STACK *stack; + QUAD_EDGE edge; +{ + QUAD_EDGE nex; + nex = onext(edge); + while (nex != edge) { + if (seen(nex) == 0) { + seen(nex) = 1; + push_edge(stack, nex); + } + nex = onext(nex); + } +} + +void +push_nonzero_ring(stack, edge) + struct EDGE_STACK *stack; + QUAD_EDGE edge; +{ + QUAD_EDGE nex; + nex = onext(edge); + while (nex != edge) { + if (seen(nex) != 0) { + seen(nex) = 0; + push_edge(stack, nex); + } + nex = onext(nex); + } +} + +void +zero_seen(my_stack,edge) + QUAD_EDGE edge; + struct EDGE_STACK *my_stack; +{ + my_stack->ptr = 0; + push_nonzero_ring(my_stack, edge); + while (my_stack->ptr != 0) { + edge = pop_edge(my_stack); + push_nonzero_ring(my_stack, sym(edge)); + } +} + |