diff options
Diffstat (limited to 'test/ccured_olden/voronoi')
| -rw-r--r-- | test/ccured_olden/voronoi/.cvsignore | 14 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/CVS/Entries | 8 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/CVS/Repository | 1 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/CVS/Root | 1 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/Makefile | 72 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/README | 22 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/defines.h | 186 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/newvor.c | 711 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/output.c | 173 | ||||
| -rw-r--r-- | test/ccured_olden/voronoi/vector.c | 80 |
10 files changed, 1268 insertions, 0 deletions
diff --git a/test/ccured_olden/voronoi/.cvsignore b/test/ccured_olden/voronoi/.cvsignore new file mode 100644 index 00000000..3923b56f --- /dev/null +++ b/test/ccured_olden/voronoi/.cvsignore @@ -0,0 +1,14 @@ +*.i +*_all.c +*cil.c +*box.c +code +data.in +data.out +allcfiles +ope.m +ope.m +*_comb.c +*cured.c +*.optim.c +changes diff --git a/test/ccured_olden/voronoi/CVS/Entries b/test/ccured_olden/voronoi/CVS/Entries new file mode 100644 index 00000000..0d7c09d8 --- /dev/null +++ b/test/ccured_olden/voronoi/CVS/Entries @@ -0,0 +1,8 @@ +/Makefile/1.3/Sat Jul 21 16:09:41 2001// +/README/1.1/Tue Jun 12 04:17:59 2001// +/defines.h/1.2/Mon Jul 9 21:12:42 2001// +/newvor.c/1.2/Mon Jul 9 21:12:42 2001// +/output.c/1.2/Mon Jul 9 21:12:42 2001// +/vector.c/1.2/Mon Jul 9 21:12:42 2001// +/.cvsignore/1.5/Thu Nov 8 23:35:28 2001// +D diff --git a/test/ccured_olden/voronoi/CVS/Repository b/test/ccured_olden/voronoi/CVS/Repository new file mode 100644 index 00000000..8d79f068 --- /dev/null +++ b/test/ccured_olden/voronoi/CVS/Repository @@ -0,0 +1 @@ +cil/test/olden/voronoi diff --git a/test/ccured_olden/voronoi/CVS/Root b/test/ccured_olden/voronoi/CVS/Root new file mode 100644 index 00000000..35f411e9 --- /dev/null +++ b/test/ccured_olden/voronoi/CVS/Root @@ -0,0 +1 @@ +/home/cvs-repository diff --git a/test/ccured_olden/voronoi/Makefile b/test/ccured_olden/voronoi/Makefile new file mode 100644 index 00000000..2822178d --- /dev/null +++ b/test/ccured_olden/voronoi/Makefile @@ -0,0 +1,72 @@ +# /* For copyright information, see olden_v1.0/COPYRIGHT */ + +BINARY = voronoi.exe +PROGS = newvor vector output args ssplain trusted_voronoi + +ifdef _MSVC +CC = cl +DEF = /D +CONLY = /c +OBJOUT = /Fo +EXEOUT = /Fe +OBJ = .obj + +OPTFLAGS = /Ox +LIBS = + +else + +CC = gcc -arch ppc +CCOMP=../../../../ccomp +CCOMPFLAGS=-dump-c + +DEF = -D +CONLY = -c +OBJOUT= -o +EXEOUT= -o + +OBJ = .o + +OPTFLAGS = -g -Wall -O3 + +LIBS = +LIBPATH = +endif + +SRC = .c +ASM = .s +EXTRA_CDEFS = $(DEF)I_TIME $(DEF)I_SYS_TIME $(DEF)ULTRIX +# sm: ?? +ifdef _MSVC +CDEFS = $(DEF)PLAIN $(DEF)SS_PLAIN $(DEF)OLDEN +else +CDEFS = $(DEF)PLAIN $(DEF)OLDEN +endif +SRCS = $(addsuffix $(SRC),$(FILES)) +OBJS = $(addsuffix $(OBJ),$(FILES)) +ASMS = $(addsuffix $(ASM),$(FILES)) + +all_s: $(PROGS:%=%.s) + +all: $(PROGS:%=%.compcert) + +all_gcc: $(PROGS:%=%.gcc) + +%.compcert: %.s + $(CC) $(CFLAGS) $(LDFALGS) $(OPTFLAGS) -o $*.compcert $*.s $(LIBS) + +%.s: %.c ../../../../ccomp + $(CCOMP) $(CCOMPFLAGS) $(CDEFS) $(EXTRA_CDEFS) $(MY_CDEFS) $*.c + +%.gcc: %.c + $(CC) $(CFLAGS) $(LDFALGS) $(OPTFLAGS) -o $*.gcc $*.c $(LIBS) + +#$(BINARY): $(OBJS) +# $(CC) $(LDFALGS) $(OPTFLAGS) $(EXEOUT)$@ $(OBJS) $(LIBPATH) $(LIBS) + +#%$(OBJ) : %$(SRC) +# $(CC) $(CDEFS) $(EXTRA_CDEFS) $(MY_CDEFS) $(OPTFLAGS) $(CONLY) $< + +clean: + rm -f $(BINARY) $(OBJS) *~ + diff --git a/test/ccured_olden/voronoi/README b/test/ccured_olden/voronoi/README new file mode 100644 index 00000000..701a0134 --- /dev/null +++ b/test/ccured_olden/voronoi/README @@ -0,0 +1,22 @@ +/* For copyright information, see olden_v1.01/COPYRIGHT */ +********************** +olden_v1.01/benchmarks/voronoi/README +June 1996 +Martin C. Carlisle + +this directory contains the Voronoi Diagram benchmark: + +L. Guibas and J. Stolfi. "General Subdivisions and Voronoi Diagrams" +ACM Trans. on Graphics 4(2):74-123, 1985. + +Adapted for Olden by Martin C. Carlisle + +********************** + +Makefile - use "make voronoi" to create executable + +args.c - process command line args +vector.c - vector math stuff +newvor.c - main routine +defines.h - definitions +output.c - used for outputting results. diff --git a/test/ccured_olden/voronoi/defines.h b/test/ccured_olden/voronoi/defines.h new file mode 100644 index 00000000..651336db --- /dev/null +++ b/test/ccured_olden/voronoi/defines.h @@ -0,0 +1,186 @@ +/* For copyright information, see olden_v1.0/COPYRIGHT */ + +extern double sqrt(); +extern double exp(); +extern double log(); +extern double drand(); + +#include <stdio.h> +#include <stdlib.h> + +#ifndef SS_PLAIN +#include "mem-ref.h" +#include "future-cell.h" + +#define NULL 0 +#endif SS_PLAIN + +#define NEW 1 +#define EPSILON (1.0e-6) + +#define BOOLEAN int +#ifndef TRUE +#define TRUE 1 +#endif TRUE +#define FALSE 0 + +struct edge_rec +{ + struct VERTEX *v; + struct edge_rec *next; + int wasseen; + int more_data; /* 16 byte align this thing */ +}; + +struct get_point +{ + struct VERTEX *v; + double curmax; + int seed; +}; + +#ifdef NEW +typedef struct edge_rec *EDGE_PTR; +typedef struct VERTEX *VERTEX_PTR; +typedef EDGE_PTR QUAD_EDGE; +#else +typedef int EDGE_PTR; +typedef int VERTEX_PTR; +#endif + +struct VEC2 { + double x,y; + double norm; +}; + +struct VERTEX { + struct VEC2 v; + struct VERTEX *left, *right; +}; + +typedef struct +{ + QUAD_EDGE left, right; +} EDGE_PAIR; + +#ifdef NEW +#define onext(a) (a)->next +#else +#define onext(a) next[a] +#endif +#define oprev(a) rot(onext(rot(a))) +#define lnext(a) rot(onext(rotinv(a))) +#define lprev(a) sym(onext(a)) +#define rnext(a) rotinv(onext(rot(a))) +#define rprev(a) onext(sym(a)) +#define dnext(a) sym(onext(sym(a))) +#define dprev(a) rotinv(onext(rotinv(a))) + +#ifdef NEW +#define X(r) r->v.x +#define Y(r) r->v.y +#define NORM(r) r->v.norm +#else +#define X(a) va[a].v.x +#define Y(a) va[a].v.y +#define NORM(a) va[a].norm +#endif +#ifdef NEW +#define orig(a) (a)->v +#define dest(a) orig(sym(a)) +#define seen(a) (a)->wasseen +#else +#define orig(a) org[a] +#define dest(a) orig(sym(a)) +#define seen(a) see[a] +#endif + +#ifndef NEW +#define origv(a) va[orig(a)].v +#define destv(a) va[dest(a)].v +#else +#define origv(a) orig(a)->v +#define destv(a) dest(a)->v +#endif + +#define ALLOC_SIZE (sizeof(struct edge_rec)) +#ifndef PLAIN +#define SIZE (sizeof(struct edge_rec) << PN_BITS) +#define ANDF ((4*sizeof(struct edge_rec) - 1) << PN_BITS) +#else +#define SIZE sizeof(struct edge_rec) +#define ANDF (4*sizeof(struct edge_rec) - 1) +#endif + +#define sym(a) ((QUAD_EDGE) (((unsigned int) (a)) ^ 2*SIZE)) +#define rot(a) ((QUAD_EDGE) ( (((unsigned int) (a) + 1*SIZE) & ANDF) | ((unsigned int) (a) & ~ANDF) )) +#define rotinv(a) ((QUAD_EDGE) ( (((unsigned int) (a) + 3*SIZE) & ANDF) | ((unsigned int) (a) & ~ANDF) )) +#define base(a) ((QUAD_EDGE) ((unsigned int a) & ~ANDF)) + +struct EDGE_STACK { + int ptr; + QUAD_EDGE *elts ; + int stack_size; +}; + +extern QUAD_EDGE alloc_edge(); +extern QUAD_EDGE makeedge(); +extern void splice(); +extern void swapedge(); +extern void deleteedge(); +extern QUAD_EDGE build_delaunay_triangulation(); +extern EDGE_PAIR build_delaunay(); +extern EDGE_PAIR do_merge(); +extern QUAD_EDGE connect_left(); +extern QUAD_EDGE connect_right(); + +extern void zero_seen(); +extern QUAD_EDGE pop_edge(); + +#ifdef SS_PLAIN +#define drand(seed) (((double) (seed=olden_random(seed))) / (double) 2147483647) +#else SS_PLAIN +#define drand(seed) (((double) (seed=random(seed))) / (double) 2147483647) +#endif SS_PLAIN + +#ifdef DEFINE_GLOBALS +#define EXTERN +#else +#define EXTERN extern +#endif + +EXTERN VERTEX_PTR *vp ; +EXTERN struct VERTEX *va ; +EXTERN EDGE_PTR *next ; +EXTERN VERTEX_PTR *org ; +EXTERN int num_vertices, num_edgeparts, stack_size ; +EXTERN int to_lincoln , to_off, to_3d_out, to_color , voronoi , delaunay , interactive , ahost ; +EXTERN char *see; + +#ifdef ODEFINE_GLOBALS +#define OEXTERN +#else +#define OEXTERN extern +#endif + +#define my_alloc(str_name, str_type, str_cnt) \ + if (NULL == (str_name = (str_type *) \ + mymalloc((unsigned) (str_cnt ) * (unsigned) (sizeof(str_type)))))\ + exit(printf(" cannot malloc (str_name) \n")) + +#define VERTEX_ALLOC 1000 + +#define RED 123 +#define GREEN (RED + 1) + +#define CY_SOLID 1 +#define CY_DOTTED 2 + + +int ccw(VERTEX_PTR a , VERTEX_PTR b , VERTEX_PTR c ); +int olden_random(int seed ); +void filestuff(void); +int dealwithargs(int argc , char * argv[] ); +void output_voronoi_diagram(QUAD_EDGE edge , int nv , struct EDGE_STACK * + my_stack ); + 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)); + } +} + diff --git a/test/ccured_olden/voronoi/output.c b/test/ccured_olden/voronoi/output.c new file mode 100644 index 00000000..a9331f99 --- /dev/null +++ b/test/ccured_olden/voronoi/output.c @@ -0,0 +1,173 @@ +/* For copyright information, see olden_v1.0/COPYRIGHT */ + +#define ODEFINE_GLOBALS +#include "defines.h" + +#ifdef SS_PLAIN +#include "ssplain.h" +#endif SS_PLAIN + +extern struct VEC2 V2_sum(); +extern struct VEC2 V2_sub(); +extern struct VEC2 V2_times(); +extern double V2_cprod(); +extern struct VEC2 V2_cross(); +extern double V2_dot(); +extern double V2_magn(); + +/****************************************************************/ +/* Voronoi Output Routines */ +/****************************************************************/ + +void +plot_dedge(p1, p2) +VERTEX_PTR p1, p2; +{ + double x1,x2,y1,y2; + + x1=X(p1); + y1=Y(p1); + x2=X(p2); + y2=Y(p2); + /* plots a Delaunay-triangulation edge on your favorite device. */ + chatting("Dedge %g %g %g %g \n",(float) x1, (float) y1, + (float) x2, (float) y2); +} + +void +plot_vedge(p1, p2) + struct VEC2 p1, p2; +{ + /* plots a Voronoi-diagram edge on your favorite device. */ + chatting("Vedge %g %g %g %g \n",(float) p1.x, (float) p1.y, (float) p2.x, + (float) p2.y); +} + +struct VEC2 circle_center(a,b,c) + /* returns the center of the circle passing through A, B & C. */ + struct VEC2 a,b,c; +{ + struct VEC2 p; + double d1,d2,d3,d4; + struct VEC2 vv1, vv2, vv3, vv4, vv5, vv6, vv7 ; + vv1 = V2_sub(b,c); + d1 = V2_magn( vv1 ); + vv1 = V2_sum(a,b); + vv2 = V2_times(0.5, vv1); + if (d1 < 0.0) /*there is no intersection point, the bisectors coincide. */ + return(vv2); + else { + vv3 = V2_sub(b,a); + vv4 = V2_sub(c,a); + d3 = V2_cprod(vv3, vv4) ; + d2 = -2.0 * d3 ; + vv5 = V2_sub(c,b); + d4 = V2_dot(vv5, vv4); + vv6 = V2_cross(vv3); + vv7 = V2_times(d4/ d2 , vv6); + p = V2_sum(vv2, vv7); + return(p); + } +} + +int *earray; + +void +output_voronoi_diagram(edge,nv,my_stack) + QUAD_EDGE edge; + int nv; + struct EDGE_STACK *my_stack; +{ + QUAD_EDGE nex, prev, snex, sprev; + struct VEC2 cvxvec, center; + double ln; + double d1; + struct VEC2 vv1, vv2, vv3; + + if (voronoi) { + zero_seen(my_stack,edge); + nex = edge; + /* Plot VD edges with one endpoint at infinity. */ + do { + struct VEC2 v21,v22,v23; + QUAD_EDGE tmp; + + v21=destv(nex); + v22=origv(nex); + tmp=onext(nex); + v23=destv(tmp); + cvxvec = V2_sub(v21,v22/*destv(nex), origv(nex)*/); + center = circle_center(v22,v21,v23 + /*origv(nex), destv(nex), destv(onext(nex))*/); + vv1 = V2_sum(v22,v21/*origv(nex), destv(nex)*/) ; + vv2 = V2_times(0.5, vv1); + vv3 = V2_sub(center, vv2) ; + ln = 1.0 + V2_magn(vv3); + d1 = ln/V2_magn(cvxvec); + vv1 = V2_cross(cvxvec); + vv2 = V2_times(d1, vv1) ; + vv3 = V2_sum(center, vv2); + plot_vedge( center, vv3 ) ; + nex = rnext(nex); + } while (nex != edge); + } + + /* Plot DT edges and finite VD edges. */ + + my_stack->ptr = 0; + push_ring(my_stack, edge); + chatting("mystack_ptr = %d\n",my_stack->ptr); + while (my_stack->ptr != 0) { + VERTEX_PTR v1,v2,v3,v4; + double d1,d2; + + + edge = pop_edge(my_stack); + if (seen(edge) == 1) + { + { + prev = edge; + nex = onext(edge); + do { + v1=orig(prev); + d1=X(v1); + v2=dest(prev); + d2=X(v2); + if (d1 >= d2 /*X(orig(prev)) >= X(dest(prev))*/) + { + /*plot_dedge(orig(prev), dest(prev));*/ + plot_dedge(v1,v2); + sprev = sym(prev); + snex = onext(sprev); + v1=orig(prev); + v2=dest(prev); + v3=dest(nex); + v4=dest(snex); + if (ccw(v1,v2,v3/*orig(prev), dest(prev), dest(nex)*/) + != ccw(v1,v2,v4/*orig(prev),dest(prev),dest(snex)*/)) + { + struct VEC2 v21, v22, v23; + v21 = origv(prev); + v22 = destv(prev); + v23 = destv(nex); + vv1 = circle_center(v21,v22,v23/*origv(prev), + destv(prev),destv(nex)*/); + v21 = origv(sprev); + v22 = destv(sprev); + v23 = destv(snex); + vv2 = circle_center(v21,v22,v23/*origv(sprev), + destv(sprev),destv(snex)*/); + plot_vedge( vv1 , vv2 ); + } + } + seen(prev) = 2; + prev = nex; + nex = onext(nex); + } while (prev != edge); + } + } + push_ring(my_stack, sym(edge)); + } + zero_seen(my_stack, edge); +} + diff --git a/test/ccured_olden/voronoi/vector.c b/test/ccured_olden/voronoi/vector.c new file mode 100644 index 00000000..1a288d7b --- /dev/null +++ b/test/ccured_olden/voronoi/vector.c @@ -0,0 +1,80 @@ +/* For copyright information, see olden_v1.0/COPYRIGHT */ + +#include "defines.h" + +/****************************************************************/ +/* Vector Routines. */ +/* From CMU vision library. */ +/* They are used only for the VD, not the DT. */ +/* They are slow because of large call-by-value parameters.*/ +/****************************************************************/ + +/* V2_cprod: forms triple scalar product of [u,v,k], where k = u cross v */ +/* (returns the magnitude of u cross v in space)*/ + +double V2_cprod(u,v) + struct VEC2 u,v; +{ + return(u.x * v.y - u.y * v.x); +} + + +/* V2_dot: vector dot product */ + +double V2_dot(u,v) +struct VEC2 u,v; +{ + return(u.x * v.x + u.y * v.y); +} + +/* V2_times: multiply a vector by a scalar */ + +struct VEC2 V2_times(c,v) + double c; + struct VEC2 v; +{ + struct VEC2 ans; + ans.x = c * v.x; + ans.y = c * v.y; + return(ans); +} + +/* V2_sum, V2_sub: Vector addition and subtraction */ + +struct VEC2 V2_sum(u,v) + struct VEC2 u,v; +{ + struct VEC2 ans; + + ans.x = u.x + v.x; + ans.y = u.y + v.y; + return(ans); +} + +struct VEC2 V2_sub(u,v) + struct VEC2 u,v; +{ + struct VEC2 ans; + ans.x = u.x - v.x; + ans.y = u.y - v.y; + return(ans); +} + +/* V2_magn: magnitude of vector */ + +double V2_magn(u) + struct VEC2 u; +{ + return(sqrt(u.x*u.x+u.y*u.y)); +} + +/* returns k X v (cross product). this is a vector perpendicular to v */ + +struct VEC2 V2_cross(v) + struct VEC2 v; +{ + struct VEC2 ans; + ans.x = v.y; + ans.y = -v.x; + return(ans); +} |