From 285f5bec5bb03d4e825e5d866e94008088dd6155 Mon Sep 17 00:00:00 2001 From: xleroy Date: Sat, 9 Aug 2008 08:06:33 +0000 Subject: Ajout nouveaux tests git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@708 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e --- test/raytracer/.depend | 20 ++ test/raytracer/Makefile | 43 +++ test/raytracer/arrays.c | 30 ++ test/raytracer/arrays.h | 31 +++ test/raytracer/config.h | 27 ++ test/raytracer/eval.c | 672 +++++++++++++++++++++++++++++++++++++++++++++ test/raytracer/eval.h | 18 ++ test/raytracer/gml.h | 73 +++++ test/raytracer/gmllexer.c | 297 ++++++++++++++++++++ test/raytracer/gmllexer.h | 21 ++ test/raytracer/gmlparser.c | 102 +++++++ test/raytracer/gmlparser.h | 3 + test/raytracer/intersect.c | 416 ++++++++++++++++++++++++++++ test/raytracer/intersect.h | 11 + test/raytracer/kal.gml | 78 ++++++ test/raytracer/light.c | 126 +++++++++ test/raytracer/light.h | 33 +++ test/raytracer/main.c | 13 + test/raytracer/matrix.c | 102 +++++++ test/raytracer/matrix.h | 18 ++ test/raytracer/memory.c | 60 ++++ test/raytracer/object.c | 214 +++++++++++++++ test/raytracer/object.h | 36 +++ test/raytracer/point.h | 18 ++ test/raytracer/render.c | 128 +++++++++ test/raytracer/render.h | 9 + test/raytracer/simplify.c | 177 ++++++++++++ test/raytracer/simplify.h | 3 + test/raytracer/surface.c | 148 ++++++++++ test/raytracer/surface.h | 7 + test/raytracer/vector.c | 74 +++++ test/raytracer/vector.h | 23 ++ 32 files changed, 3031 insertions(+) create mode 100644 test/raytracer/.depend create mode 100644 test/raytracer/Makefile create mode 100644 test/raytracer/arrays.c create mode 100644 test/raytracer/arrays.h create mode 100644 test/raytracer/config.h create mode 100644 test/raytracer/eval.c create mode 100644 test/raytracer/eval.h create mode 100644 test/raytracer/gml.h create mode 100644 test/raytracer/gmllexer.c create mode 100644 test/raytracer/gmllexer.h create mode 100644 test/raytracer/gmlparser.c create mode 100644 test/raytracer/gmlparser.h create mode 100644 test/raytracer/intersect.c create mode 100644 test/raytracer/intersect.h create mode 100644 test/raytracer/kal.gml create mode 100644 test/raytracer/light.c create mode 100644 test/raytracer/light.h create mode 100644 test/raytracer/main.c create mode 100644 test/raytracer/matrix.c create mode 100644 test/raytracer/matrix.h create mode 100644 test/raytracer/memory.c create mode 100644 test/raytracer/object.c create mode 100644 test/raytracer/object.h create mode 100644 test/raytracer/point.h create mode 100644 test/raytracer/render.c create mode 100644 test/raytracer/render.h create mode 100644 test/raytracer/simplify.c create mode 100644 test/raytracer/simplify.h create mode 100644 test/raytracer/surface.c create mode 100644 test/raytracer/surface.h create mode 100644 test/raytracer/vector.c create mode 100644 test/raytracer/vector.h (limited to 'test/raytracer') diff --git a/test/raytracer/.depend b/test/raytracer/.depend new file mode 100644 index 0000000..c31fbcd --- /dev/null +++ b/test/raytracer/.depend @@ -0,0 +1,20 @@ +arrays.o: arrays.c config.h arrays.h +eval.o: eval.c config.h arrays.h gml.h point.h vector.h eval.h object.h \ + light.h render.h +gmllexer.o: gmllexer.c config.h gmllexer.h gml.h +gmlparser.o: gmlparser.c config.h arrays.h gml.h gmllexer.h gmlparser.h +intersect.o: intersect.c config.h arrays.h point.h vector.h eval.h \ + object.h matrix.h intersect.h +light.o: light.c config.h point.h vector.h eval.h object.h intersect.h \ + light.h +main.o: main.c config.h arrays.h gmllexer.h gmlparser.h eval.h +matrix.o: matrix.c config.h point.h vector.h matrix.h +memory.o: memory.c config.h +object.o: object.c config.h arrays.h point.h vector.h eval.h object.h \ + matrix.h +render.o: render.c config.h point.h vector.h matrix.h eval.h object.h \ + light.h intersect.h surface.h simplify.h render.h +simplify.o: simplify.c config.h point.h vector.h arrays.h matrix.h eval.h \ + object.h simplify.h +surface.o: surface.c config.h point.h vector.h eval.h object.h surface.h +vector.o: vector.c config.h point.h vector.h diff --git a/test/raytracer/Makefile b/test/raytracer/Makefile new file mode 100644 index 0000000..cd442c5 --- /dev/null +++ b/test/raytracer/Makefile @@ -0,0 +1,43 @@ +CC=../../ccomp +CFLAGS=-U__GNUC__ -stdlib ../../runtime -dclight -dasm +LIBS= + +OBJS=memory.o gmllexer.o gmlparser.o eval.o \ + arrays.o vector.o matrix.o object.o intersect.o surface.o light.o \ + simplify.o render.o main.o + +render: $(OBJS) + $(CC) $(CFLAGS) -o render $(OBJS) $(LIBS) + +clean: + rm -f *.o *.light.c *.s *.ppm render + +include .depend + +depend: + gcc -MM *.c > .depend + +gcc0: + $(MAKE) clean + $(MAKE) CC=gcc CFLAGS="-O0 -Wall" + mv render render.gcc0 + $(MAKE) clean + +gcc1: + $(MAKE) clean + $(MAKE) CC=gcc CFLAGS="-O1 -Wall" + mv render render.gcc1 + $(MAKE) clean + +gcc2: + $(MAKE) clean + $(MAKE) CC=gcc CFLAGS="-O2 -Wall" + mv render render.gcc2 + $(MAKE) clean + +gcc3: + $(MAKE) clean + $(MAKE) CC=gcc CFLAGS="-O3 -Wall" + mv render render.gcc3 + $(MAKE) clean + diff --git a/test/raytracer/arrays.c b/test/raytracer/arrays.c new file mode 100644 index 0000000..43508ee --- /dev/null +++ b/test/raytracer/arrays.c @@ -0,0 +1,30 @@ +#include "config.h" +#include "arrays.h" + +struct array * alloc_array(int eltsize, int initialsize) +{ + struct array * a = arena_alloc(sizeof(struct array)); + a->size = 0; + a->capacity = initialsize; + a->data = arena_alloc(eltsize * initialsize); + return a; +} + +void grow_array(int eltsize, struct array * arr) +{ + void * newdata = arena_alloc(arr->capacity * 2 * eltsize); + memcpy(newdata, arr->data, arr->size * eltsize); + arr->data = newdata; + arr->capacity *= 2; +} + +struct array * copy_array(int eltsize, struct array * arr, int extrasize) +{ + struct array * a = arena_alloc(sizeof(struct array)); + a->size = arr->size; + a->capacity = arr->size + extrasize; + a->data = arena_alloc(eltsize * a->capacity); + memcpy(a->data, arr->data, eltsize * a->size); + return a; +} + diff --git a/test/raytracer/arrays.h b/test/raytracer/arrays.h new file mode 100644 index 0000000..47634f6 --- /dev/null +++ b/test/raytracer/arrays.h @@ -0,0 +1,31 @@ +/* Resizable arrays */ + +struct array { + int size; /* Number of elements in the array */ + int capacity; /* Max number of elements before resizing */ + void * data; /* Actual data -- variable sized */ +}; + +struct array * alloc_array(int eltsize, int initialsize); + +void grow_array(int eltsize, struct array * arr); + +struct array * copy_array(int eltsize, struct array * arr, int extrasize); + +#define new_array(ty,sz) alloc_array(sizeof(ty), sz) + +#define data_array(ty,arr) ((ty *) (arr)->data) + +#define get_array(ty,arr,idx) (data_array(ty,arr)[idx]) + +#define get_array_large(dst,ty,arr,idx) \ + ASSIGN(dst, data_array(ty,arr)[idx]) + +#define set_array(ty,arr,idx,newval) (data_array(ty,arr)[idx] = (newval)) + +#define set_array_large(ty,arr,idx,newval) \ + ASSIGN(data_array(ty,arr)[idx], newval) + +#define extend_array(ty,arr) \ + if ((arr)->size >= (arr)->capacity) grow_array(sizeof(ty), (arr)); \ + (arr)->size++ diff --git a/test/raytracer/config.h b/test/raytracer/config.h new file mode 100644 index 0000000..ed5b847 --- /dev/null +++ b/test/raytracer/config.h @@ -0,0 +1,27 @@ +#include +#include +#include +#include +#include +#include + +#ifdef SINGLE_PRECISION +typedef float flt; +#else +typedef double flt; +#endif + +#if 0 +extern void abord(void); + +#define assert(cond) \ + if (!(cond)) { fprintf(stderr, "%s:%u: failed assertion\n", __FILE__, __LINE__); abort(); } + +#endif + +#define ASSIGN(lv,rv) memcpy(&(lv), &(rv), sizeof(rv)) + +void arena_init(void); +void arena_clear(void); +void * arena_alloc(int size); + diff --git a/test/raytracer/eval.c b/test/raytracer/eval.c new file mode 100644 index 0000000..4d3f3dc --- /dev/null +++ b/test/raytracer/eval.c @@ -0,0 +1,672 @@ +/* Evaluator for GML */ + +#include "config.h" +#include "arrays.h" +#include "gml.h" +#include "point.h" +#include "vector.h" +#include "eval.h" +#include "object.h" +#include "light.h" +#include "render.h" + +struct value { + enum { B, I, R, S, Arr, Clos, Point, Obj, Light } tag; + union { + int i; /* B, I */ + flt r; /* R */ + char * s; /* S */ + struct array * arr; /* Arr */ + struct closure clos; /* Clos */ + struct point * point; /* Point */ + struct object * obj; /* Obj */ + struct light * light; /* Light */ + } u; +}; + +struct binding { + char * name; + int mutable; + struct value val; +}; + +/* Lookup an identifier in an environment */ + +static int lookup(struct array * env, char * name, /*out*/ struct value * res) +{ + int i; + for (i = env->size - 1; i >= 0; i--) { + struct binding * b = &get_array(struct binding, env, i); + if (name == b->name) { + ASSIGN(*res, b->val); + return 1; + } + } + return 0; +} + +/* Assign an identifier in an environment */ + +static void assign(struct array * env, char * name, struct value * newval) +{ + int i; + struct binding * b; + for (i = env->size - 1; i >= 0; i--) { + b = &get_array(struct binding, env, i); + if (! b->mutable) break; + if (name == b->name) { + ASSIGN(b->val, *newval); + return; + } + } + extend_array(struct binding, env); + b = &get_array(struct binding, env, env->size - 1); + b->name = name; + b->mutable = 1; + ASSIGN(b->val, *newval); +} + +/* Take an immutable copy of an environment */ + +static struct array * snapshot_env(struct array * env) +{ + int i; + struct array * nenv = copy_array(sizeof(struct binding), env, 10); + for (i = 0; i < nenv->size; i++) + get_array(struct binding, nenv, i).mutable = 0; + return nenv; +} + +/* Utility math functions */ + +static inline flt degrees_to_radians(flt d) +{ return d * (M_PI / 180.0); } + +static inline flt radians_to_degrees(flt d) +{ return d * (180.0 / M_PI); } + +static inline flt deg_cos(flt a) +{ return cos(degrees_to_radians(a)); } + +static inline flt deg_sin(flt a) +{ return sin(degrees_to_radians(a)); } + +static inline flt deg_acos(flt a) +{ return radians_to_degrees(acos(a)); } + +static inline flt deg_asin(flt a) +{ return radians_to_degrees(asin(a)); } + +static inline flt clampf(flt r) +{ + if (r < 0.0) return 0.0; + if (r > 1.0) return 1.0; + return r; +} + +/* For error printing */ + +static char * operator_names [] = { + "", "", "", "", "", + "", "", "", "acos", "addi", "addf", "apply", + "asin", "clampf", "cone", "cos", "cube", "cylinder", "difference", + "divi", "divf", "eqi", "eqf", "floor", "frac", "get", "getx", "gety", + "getz", "if", "intersect", "length", "lessi", "lessf", "light", + "modi", "muli", "mulf", "negi", "negf", "plane", "point", + "pointlight", "real", "render", "rotatex", "rotatey", "rotatez", + "scale", "sin", "sphere", "spotlight", "sqrt", "subi", "subf", + "translate", "union", "uscale", "print", +}; + +/* Convert a GML array of lights into a C array of lights */ + +static struct light ** light_array(struct array * arr) +{ + int sz = arr->size; + struct light ** l = arena_alloc((sz + 1) * sizeof(struct light *)); + int i; + for (i = 0; i < sz; i++) { + struct value * v = &get_array(struct value, arr, i); + if (v->tag != Light) { + fprintf(stderr, "Light expected in array argument to `render'\n"); + exit(2); + } + l[i] = v->u.light; + } + l[sz] = NULL; + return l; +} + +/* Pretty-print a value */ + +static void print_value(struct value * s) +{ + switch (s->tag) { + case B: + printf("%s\n", s->u.i ? "true" : "false"); break; + case I: + printf("%d\n", s->u.i); break; + case R: + printf("%e\n", s->u.r); break; + case S: + printf("\"%s\"\n", s->u.s); break; + case Arr: + printf("array\n"); break; + case Clos: + printf("closure\n"); break; + case Point: + printf("point %e %e %e\n", + s->u.point->x, s->u.point->y, s->u.point->z); + break; + case Obj: + printf("object\n"); break; + case Light: + printf("light\n"); break; + } +} + +/* The evaluation stack */ + +#define MAIN_STACK_SIZE 10000 +#define SURFACE_STACK_SIZE 1000 + +static struct value main_stack[MAIN_STACK_SIZE]; +static struct value surface_stack[SURFACE_STACK_SIZE]; + +/* Error handling functions */ +#define ERRORFUN(name, msg) \ + static void name(void) { fprintf(stderr, msg); exit(2); } + +ERRORFUN(stack_overflow, "Stack overflow\n") +ERRORFUN(stack_underflow, "Stack underflow\n") +ERRORFUN(type_error, "Type error\n") +ERRORFUN(division_by_zero, "Division by zero\n") +ERRORFUN(bound_error, "Out-of-bound array access\n") +ERRORFUN(negative_sqrt, "Square root of negative number\n") + +/* Macros for stack checking and type checking */ + +#define push() if (--sp < bos) stack_overflow() +#define can_pop(n) if (sp + (n) > tos) stack_underflow() +#define check(n,ty) if (sp[n].tag != ty) type_error() + +/* Execute the given token list in the given environment */ + +static struct value * execute_list(struct array * code, + struct array * env, + struct value * bos, + struct value * tos, + struct value * sp) +{ + int i; + struct tok * t; + + for (i = 0; i < code->size; i++) { + t = &get_array(struct tok, code, i); + switch (t->tag) { + case Identifier: + push(); + if (! lookup(env, t->u.s, sp)) { + fprintf(stderr, "Unbound identifier %s\n", t->u.s); + exit(2); + } + break; + case Binder: + can_pop(1); + assign(env, t->u.s, sp); + sp++; + break; + case Boolean: + push(); + sp[0].tag = B; sp[0].u.i = t->u.i; + break; + case Integer: + push(); + sp[0].tag = I; sp[0].u.i = t->u.i; + break; + case Real: + push(); + sp[0].tag = R; sp[0].u.r = t->u.d; + break; + case String: + push(); + sp[0].tag = S; sp[0].u.s = t->u.s; + break; + case Array: { + struct value * esp = execute_list(t->u.a, env, bos, sp, sp); + int sz = sp - esp; + struct array * a = new_array(struct value, sz); + int j; + a->size = sz; + for (j = 0; j < sz; j++) set_array_large(struct value, a, j, sp[-1-j]); + push(); + sp[0].tag = Arr; sp[0].u.arr = a; + break; + } + case Function: + push(); + sp[0].tag = Clos; + sp[0].u.clos.code = t->u.a; + sp[0].u.clos.env = snapshot_env(env); + break; + case Op_acos: + can_pop(1); + check(0, R); + sp[0].u.r = deg_acos(sp[0].u.r); + break; + case Op_addi: + can_pop(2); + check(1, I); + check(0, I); + sp[1].u.i = sp[1].u.i + sp[0].u.i; + sp += 1; + break; + case Op_addf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].u.r = sp[1].u.r + sp[0].u.r; + sp += 1; + break; + case Op_apply: + can_pop(1); + check(0, Clos); + sp = execute_list(sp[0].u.clos.code, sp[0].u.clos.env, bos, tos, sp + 1); + break; + case Op_asin: + can_pop(1); + check(0, R); + sp[0].u.r = deg_asin(sp[0].u.r); + break; + case Op_clampf: + can_pop(1); + check(0, R); + sp[0].u.r = clampf(sp[0].u.r); + break; + case Op_cone: + can_pop(1); + check(0, Clos); + sp[0].tag = Obj; + sp[0].u.obj = cone(&sp[0].u.clos); + break; + case Op_cos: + can_pop(1); + check(0, R); + sp[0].u.r = deg_cos(sp[0].u.r); + break; + case Op_cube: + can_pop(1); + check(0, Clos); + sp[0].tag = Obj; + sp[0].u.obj = cube(&sp[0].u.clos); + break; + case Op_cylinder: + can_pop(1); + check(0, Clos); + sp[0].tag = Obj; + sp[0].u.obj = cylinder(&sp[0].u.clos); + break; + case Op_difference: + can_pop(2); + check(1, Obj); + check(0, Obj); + sp[1].u.obj = odifference(sp[1].u.obj, sp[0].u.obj); + sp += 1; + break; + case Op_divi: + can_pop(2); + check(1, I); + check(0, I); + if (sp[0].u.i == 0) division_by_zero(); + sp[1].u.i = sp[1].u.i / sp[0].u.i; + sp += 1; + break; + case Op_divf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].u.r = sp[1].u.r / sp[0].u.r; + sp += 1; + break; + case Op_eqi: + can_pop(2); + check(1, I); + check(0, I); + sp[1].tag = B; + sp[1].u.i = (sp[1].u.i == sp[0].u.i); + sp += 1; + break; + case Op_eqf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].tag = B; + sp[1].u.i = (sp[1].u.r == sp[0].u.r); + sp += 1; + break; + case Op_floor: + can_pop(1); + check(0, R); + sp[0].tag = I; + sp[0].u.i = (int) floor(sp[0].u.r); + break; + case Op_frac: { + double rem; + can_pop(1); + check(0, R); + sp[0].u.r = modf(sp[0].u.r, &rem); + break; + } + case Op_get: { + struct array * a; + int idx; + can_pop(2); + check(1, Arr); + check(0, I); + a = sp[1].u.arr; + idx = sp[0].u.i; + if (idx < 0 || idx >= a->size) bound_error(); + get_array_large(sp[1], struct value, a, idx); + sp++; + break; + } + case Op_getx: + can_pop(1); + check(0, Point); + sp[0].tag = R; + sp[0].u.r = sp[0].u.point->x; + break; + case Op_gety: + can_pop(1); + check(0, Point); + sp[0].tag = R; + sp[0].u.r = sp[0].u.point->y; + break; + case Op_getz: + can_pop(1); + check(0, Point); + sp[0].tag = R; + sp[0].u.r = sp[0].u.point->z; + break; + case Op_if: + can_pop(3); + check(2, B); + check(1, Clos); + check(0, Clos); + if (sp[2].u.i) + sp = execute_list(sp[1].u.clos.code, sp[1].u.clos.env, bos, tos, sp + 3); + else + sp = execute_list(sp[0].u.clos.code, sp[0].u.clos.env, bos, tos, sp + 3); + break; + case Op_intersect: + can_pop(2); + check(1, Obj); + check(0, Obj); + sp[1].u.obj = ointersect(sp[1].u.obj, sp[0].u.obj); + sp += 1; + break; + case Op_length: + can_pop(1); + check(0, Arr); + sp[0].tag = I; + sp[0].u.i = sp[0].u.arr->size; + break; + case Op_lessi: + can_pop(2); + check(1, I); + check(0, I); + sp[1].tag = B; + sp[1].u.i = (sp[1].u.i < sp[0].u.i); + sp += 1; + break; + case Op_lessf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].tag = B; + sp[1].u.i = (sp[1].u.r < sp[0].u.r); + sp += 1; + break; + case Op_light: + can_pop(2); + check(1, Point); + check(0, Point); + sp[1].tag = Light; + sp[1].u.light = dirlight(sp[1].u.point, sp[0].u.point); + sp += 1; + break; + case Op_modi: + can_pop(2); + check(1, I); + check(0, I); + if (sp[0].u.i == 0) division_by_zero(); + sp[1].u.i = sp[1].u.i % sp[0].u.i; + sp += 1; + break; + case Op_muli: + can_pop(2); + check(1, I); + check(0, I); + sp[1].u.i = sp[1].u.i * sp[0].u.i; + sp += 1; + break; + case Op_mulf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].u.r = sp[1].u.r * sp[0].u.r; + sp += 1; + break; + case Op_negi: + can_pop(1); + check(0, I); + sp[0].u.i = - sp[0].u.i; + break; + case Op_negf: + can_pop(1); + check(0, R); + sp[0].u.r = - sp[0].u.r; + break; + case Op_plane: + can_pop(1); + check(0, Clos); + sp[0].tag = Obj; + sp[0].u.obj = plane(&sp[0].u.clos); + break; + case Op_point: { + struct point * p; + can_pop(3); + check(2, R); + check(1, R); + check(0, R); + p = arena_alloc(sizeof(struct point)); + p->x = sp[2].u.r; + p->y = sp[1].u.r; + p->z = sp[0].u.r; + sp += 2; + sp[0].tag = Point; + sp[0].u.point = p; + break; + } + case Op_pointlight: + can_pop(2); + check(1, Point); + check(0, Point); + sp[1].tag = Light; + sp[1].u.light = pointlight(sp[1].u.point, sp[0].u.point); + sp += 1; + break; + case Op_real: + can_pop(1); + check(0, I); + sp[0].tag = R; + sp[0].u.r = (flt) sp[0].u.i; + break; + case Op_render: + can_pop(8); + check(0, S); + check(1, I); + check(2, I); + check(3, R); + check(4, I); + check(5, Obj); + check(6, Arr); + check(7, Point); + render(sp[7].u.point, + sp[6].u.arr->size, + light_array(sp[6].u.arr), + sp[5].u.obj, + sp[4].u.i, + degrees_to_radians(sp[3].u.r), + sp[2].u.i, + sp[1].u.i, + sp[0].u.s); + sp += 8; + break; + case Op_rotatex: + can_pop(2); + check(1, Obj); + check(0, R); + sp[1].u.obj = orotatex(sp[1].u.obj, degrees_to_radians(sp[0].u.r)); + sp += 1; + break; + case Op_rotatey: + can_pop(2); + check(1, Obj); + check(0, R); + sp[1].u.obj = orotatey(sp[1].u.obj, degrees_to_radians(sp[0].u.r)); + sp += 1; + break; + case Op_rotatez: + can_pop(2); + check(1, Obj); + check(0, R); + sp[1].u.obj = orotatez(sp[1].u.obj, degrees_to_radians(sp[0].u.r)); + sp += 1; + break; + case Op_scale: + can_pop(4); + check(3, Obj); + check(2, R); + check(1, R); + check(0, R); + sp[3].u.obj = oscale(sp[3].u.obj, sp[2].u.r, sp[1].u.r, sp[0].u.r); + sp += 3; + break; + case Op_sin: + can_pop(1); + check(0, R); + sp[0].u.r = deg_sin(sp[0].u.r); + break; + case Op_sphere: + can_pop(1); + check(0, Clos); + sp[0].tag = Obj; + sp[0].u.obj = sphere(&sp[0].u.clos); + break; + case Op_spotlight: + can_pop(5); + check(4, Point); + check(3, Point); + check(2, Point); + check(1, R); + check(0, R); + sp[4].tag = Light; + sp[4].u.light = spotlight(sp[4].u.point, sp[3].u.point, sp[2].u.point, + degrees_to_radians(sp[1].u.r), sp[0].u.r); + sp += 4; + break; + case Op_sqrt: + can_pop(1); + check(0, R); + if (sp[0].u.r < 0) negative_sqrt(); + sp[0].u.r = sqrt(sp[0].u.r); + break; + case Op_subi: + can_pop(2); + check(1, I); + check(0, I); + sp[1].u.i = sp[1].u.i - sp[0].u.i; + sp += 1; + break; + case Op_subf: + can_pop(2); + check(1, R); + check(0, R); + sp[1].u.r = sp[1].u.r - sp[0].u.r; + sp += 1; + break; + case Op_translate: + can_pop(4); + check(3, Obj); + check(2, R); + check(1, R); + check(0, R); + sp[3].u.obj = otranslate(sp[3].u.obj, sp[2].u.r, sp[1].u.r, sp[0].u.r); + sp += 3; + break; + case Op_union: + can_pop(2); + check(1, Obj); + check(0, Obj); + sp[1].u.obj = ounion(sp[1].u.obj, sp[0].u.obj); + sp += 1; + break; + case Op_uscale: + can_pop(2); + check(1, Obj); + check(0, R); + sp[1].u.obj = ouscale(sp[1].u.obj, sp[0].u.r); + sp += 1; + break; + case Op_print: + can_pop(1); + print_value(sp); + sp += 1; + break; + } + } + return sp; +} + +/* Evaluate a surface function */ + +void surface_function(struct closure * clos, int face, flt u, flt v, + /*out*/ struct surface_characteristics * sc) +{ + struct value * sp; + sp = surface_stack + SURFACE_STACK_SIZE - 3; + sp[2].tag = I; + sp[2].u.i = face; + sp[1].tag = R; + sp[1].u.i = u; + sp[0].tag = R; + sp[0].u.i = v; + sp = + execute_list(clos->code, clos->env, surface_stack, + surface_stack + SURFACE_STACK_SIZE, sp); + if (sp != surface_stack + SURFACE_STACK_SIZE - 4 || + sp[0].tag != R || + sp[1].tag != R || + sp[2].tag != R || + sp[3].tag != Point) { + fprintf(stderr, "Wrong result for surface function\n"); + exit(2); + } + sc->x = sp[3].u.point->x; + sc->y = sp[3].u.point->y; + sc->z = sp[3].u.point->z; + sc->kd = sp[2].u.r; + sc->ks = sp[1].u.r; + sc->phong = sp[0].u.r; +} + +/* Execute the main program */ + +void execute_program(struct array * toklist) +{ + execute_list(toklist, new_array(struct binding, 50), + main_stack, + main_stack + MAIN_STACK_SIZE, + main_stack + MAIN_STACK_SIZE); +} diff --git a/test/raytracer/eval.h b/test/raytracer/eval.h new file mode 100644 index 0000000..610b5d9 --- /dev/null +++ b/test/raytracer/eval.h @@ -0,0 +1,18 @@ +/* Evaluator for GML */ + +struct closure { + struct array * code; + struct array * env; +}; + +struct surface_characteristics { + flt x, y, z; + flt kd; + flt ks; + flt phong; +}; + +void execute_program(struct array * toklist); + +void surface_function(struct closure * clos, int face, flt u, flt v, + /*out*/ struct surface_characteristics * sc); diff --git a/test/raytracer/gml.h b/test/raytracer/gml.h new file mode 100644 index 0000000..20ca3e9 --- /dev/null +++ b/test/raytracer/gml.h @@ -0,0 +1,73 @@ +/* The GML abstract syntax tree */ + +enum operation { + Identifier, + Binder, + Boolean, + Integer, + Real, + String, + Array, + Function, + Op_acos, + Op_addi, + Op_addf, + Op_apply, + Op_asin, + Op_clampf, + Op_cone, + Op_cos, + Op_cube, + Op_cylinder, + Op_difference, + Op_divi, + Op_divf, + Op_eqi, + Op_eqf, + Op_floor, + Op_frac, + Op_get, + Op_getx, + Op_gety, + Op_getz, + Op_if, + Op_intersect, + Op_length, + Op_lessi, + Op_lessf, + Op_light, + Op_modi, + Op_muli, + Op_mulf, + Op_negi, + Op_negf, + Op_plane, + Op_point, + Op_pointlight, + Op_real, + Op_render, + Op_rotatex, + Op_rotatey, + Op_rotatez, + Op_scale, + Op_sin, + Op_sphere, + Op_spotlight, + Op_sqrt, + Op_subi, + Op_subf, + Op_translate, + Op_union, + Op_uscale, + Op_print +}; + +struct tok { + enum operation tag; + union { + char * s; /* Identifier, Binder, String */ + int i; /* Boolean, Integer */ + flt d; /* Real */ + struct array * a; /* Array, Function */ + } u; +}; diff --git a/test/raytracer/gmllexer.c b/test/raytracer/gmllexer.c new file mode 100644 index 0000000..332126b --- /dev/null +++ b/test/raytracer/gmllexer.c @@ -0,0 +1,297 @@ +/* Lexer for GML */ + +#include "config.h" +#include "gmllexer.h" +#include "gml.h" + +#define isdigit(c) (c >= '0' && c <= '9') +#define isalpha(c) ((c >= 'A' && c <= 'Z') || (c >= 'a' && c <= 'z')) +#define isalnum(c) (isdigit(c) || isalpha(c)) +#define isprint(c) (c >= ' ' && c <= 126) + +struct lexeme current_lexeme; + +struct bucket { + int kind; + struct bucket * next; + int op; + char string[1]; +}; + +#define HASHTABLE_SIZE 256 + +static struct bucket * hashtable[HASHTABLE_SIZE] = { NULL, /*all nulls*/}; + +#define BUFFER_SIZE 1024 + +static char buffer[BUFFER_SIZE]; + +#define STORE_BUFFER(pos,c) if (pos < sizeof(buffer) - 1) buffer[pos++] = c + +static inline unsigned int hash_ident(char * s) +{ + unsigned int h; + for (h = 0; *s != 0; s++) h = (h << 4) + h + *s; + return h % HASHTABLE_SIZE; +} + +static void get_ident(int firstchar) +{ + int c, pos; + unsigned int h; + struct bucket * b; + + /* Read characters of the ident */ + buffer[0] = firstchar; + pos = 1; + while (1) { + c = getchar(); + if (! (isalnum(c) || c == '-' || c == '_')) break; + STORE_BUFFER(pos, c); + } + if (c != EOF) ungetc(c, stdin); + buffer[pos] = 0; + /* Hash the ident */ + h = hash_ident(buffer); + /* Look it up in the hash table */ + for (b = hashtable[h]; b != NULL; b = b->next) { + if (strcmp(b->string, buffer) == 0) { + /* Found: return hash table entry */ + current_lexeme.kind = b->kind; + if (b->kind == IDENTIFIER) + current_lexeme.u.s = b->string; + else + current_lexeme.u.op = b->op; + return; + } + } + /* Not found: enter it */ + b = malloc(sizeof(struct bucket) + pos); + b->kind = IDENTIFIER; + strcpy(b->string, buffer); + b->next = hashtable[h]; + hashtable[h] = b; + current_lexeme.kind = IDENTIFIER; + current_lexeme.u.s = b->string; +} + +static void get_binder(void) +{ + int c = getchar(); + if (! isalpha(c)) { + fprintf(stderr, "Bad binder /%c...\n", c); + exit(2); + } + get_ident(c); + if (current_lexeme.kind != IDENTIFIER) { + fprintf(stderr, "Binder /%s rebinds reserved identifier\n", + current_lexeme.u.s); + exit(2); + } + current_lexeme.kind = BINDER; +} + +static int get_number(int firstchar) +{ + int c, pos, is_real; + + pos = 0; + is_real = 0; + c = firstchar; + /* Initial sign */ + if (c == '-') { + STORE_BUFFER(pos, c); + c = getchar(); + if (! isdigit(c)) return -1; + } + /* Decimal number */ + do { + STORE_BUFFER(pos, c); + c = getchar(); + } while (isdigit(c)); + /* Period and fractional part */ + if (c == '.') { + is_real = 1; + STORE_BUFFER(pos, c); + c = getchar(); + if (! isdigit(c)) return -1; + do { + STORE_BUFFER(pos, c); + c = getchar(); + } while (isdigit(c)); + } + /* Exponent */ + if (c == 'e' || c == 'E') { + is_real = 1; + STORE_BUFFER(pos, c); + c = getchar(); + if (c == '-') { + STORE_BUFFER(pos, c); + c = getchar(); + } + if (! isdigit(c)) return -1; + do { + STORE_BUFFER(pos, c); + c = getchar(); + } while (isdigit(c)); + } + if (c != EOF) ungetc(c, stdin); + /* Convert to token */ + buffer[pos] = 0; + if (is_real) { + current_lexeme.kind = REAL; + current_lexeme.u.d = strtod(buffer, NULL); + } else { + current_lexeme.kind = INTEGER; + current_lexeme.u.i = atoi(buffer); + } + return 0; +} + +static int get_string() +{ + int c, pos; + pos = 0; + while (1) { + c = getchar(); + if (c == '"') break; + if (! isprint(c)) return -1; + STORE_BUFFER(pos, c); + } + buffer[pos] = 0; + current_lexeme.kind = STRING; + current_lexeme.u.s = strdup(buffer); + return 0; +} + +void get_lexeme(void) +{ + int c; + + if (current_lexeme.kind != NONE) return; + while (1) { + c = getchar(); + switch (c) { + case EOF: + current_lexeme.kind = END_OF_FILE; break; + case ' ': case '\n': case '\t': case '\r': case 11: + continue; + case '%': + do { c = getchar(); } while (c != '\n' && c != EOF); + continue; + case '/': + get_binder(); break; + case '-': case '0': case '1': case '2': case '3': case '4': + case '5': case '6': case '7': case '8': case '9': + if (get_number(c) == -1) { + fprintf(stderr, "Bad number\n"); + exit(2); + } + break; + case '"': + if (get_string() == -1) { + fprintf(stderr, "Bad string literal\n"); + exit(2); + } + break; + case '{': + current_lexeme.kind = LBRACE; break; + case '}': + current_lexeme.kind = RBRACE; break; + case '[': + current_lexeme.kind = LBRACKET; break; + case ']': + current_lexeme.kind = RBRACKET; break; + default: + if (isalpha(c)) { + get_ident(c); + } else { + fprintf(stderr, "Illegal character `%c'\n", c); + exit(2); + } + break; + } + return; + } +} + +void discard_lexeme(void) +{ + current_lexeme.kind = NONE; +} + +/* Enter the operators in the hashtable */ + +struct op_init { char * name; int kind; int op; }; + +static struct op_init operators[] = + { { "true", BOOLEAN, 1 }, + { "false", BOOLEAN, 0 }, + { "acos", OPERATOR, Op_acos }, + { "addi", OPERATOR, Op_addi }, + { "addf", OPERATOR, Op_addf }, + { "apply", OPERATOR, Op_apply }, + { "asin", OPERATOR, Op_asin }, + { "clampf", OPERATOR, Op_clampf }, + { "cone", OPERATOR, Op_cone }, + { "cos", OPERATOR, Op_cos }, + { "cube", OPERATOR, Op_cube }, + { "cylinder", OPERATOR, Op_cylinder }, + { "difference", OPERATOR, Op_difference }, + { "divi", OPERATOR, Op_divi }, + { "divf", OPERATOR, Op_divf }, + { "eqi", OPERATOR, Op_eqi }, + { "eqf", OPERATOR, Op_eqf }, + { "floor", OPERATOR, Op_floor }, + { "frac", OPERATOR, Op_frac }, + { "get", OPERATOR, Op_get }, + { "getx", OPERATOR, Op_getx }, + { "gety", OPERATOR, Op_gety }, + { "getz", OPERATOR, Op_getz }, + { "if", OPERATOR, Op_if }, + { "intersect", OPERATOR, Op_intersect }, + { "length", OPERATOR, Op_length }, + { "lessi", OPERATOR, Op_lessi }, + { "lessf", OPERATOR, Op_lessf }, + { "light", OPERATOR, Op_light }, + { "modi", OPERATOR, Op_modi }, + { "muli", OPERATOR, Op_muli }, + { "mulf", OPERATOR, Op_mulf }, + { "negi", OPERATOR, Op_negi }, + { "negf", OPERATOR, Op_negf }, + { "plane", OPERATOR, Op_plane }, + { "point", OPERATOR, Op_point }, + { "pointlight", OPERATOR, Op_pointlight }, + { "real", OPERATOR, Op_real }, + { "render", OPERATOR, Op_render }, + { "rotatex", OPERATOR, Op_rotatex }, + { "rotatey", OPERATOR, Op_rotatey }, + { "rotatez", OPERATOR, Op_rotatez }, + { "scale", OPERATOR, Op_scale }, + { "sin", OPERATOR, Op_sin }, + { "sphere", OPERATOR, Op_sphere }, + { "spotlight", OPERATOR, Op_spotlight }, + { "sqrt", OPERATOR, Op_sqrt }, + { "subi", OPERATOR, Op_subi }, + { "subf", OPERATOR, Op_subf }, + { "translate", OPERATOR, Op_translate }, + { "union", OPERATOR, Op_union }, + { "uscale", OPERATOR, Op_uscale }, + { "print", OPERATOR, Op_print }, + }; + +void init_lexer(void) +{ + int i; + for (i = 0; i < sizeof(operators) / sizeof(struct op_init); i++) { + struct op_init * opi = &(operators[i]); + struct bucket * b = malloc(sizeof(struct bucket) + strlen(opi->name)); + unsigned int h = hash_ident(opi->name); + b->kind = opi->kind; + strcpy(b->string, opi->name); + b->op = opi->op; + b->next = hashtable[h]; + hashtable[h] = b; + } +} + diff --git a/test/raytracer/gmllexer.h b/test/raytracer/gmllexer.h new file mode 100644 index 0000000..9b8d630 --- /dev/null +++ b/test/raytracer/gmllexer.h @@ -0,0 +1,21 @@ +/* Lexer for GML */ + +struct lexeme { + enum { + NONE, + OPERATOR, IDENTIFIER, BINDER, BOOLEAN, INTEGER, REAL, STRING, + LBRACE, RBRACE, LBRACKET, RBRACKET, END_OF_FILE + } kind; + union { + int op; /* for operators */ + char * s; /* for identifiers, binders, strings */ + int i; /* for integer and boolean literals */ + flt d; /* for float literals */ + } u; +}; + +extern struct lexeme current_lexeme; + +extern void get_lexeme(void); +extern void discard_lexeme(void); +extern void init_lexer(void); diff --git a/test/raytracer/gmlparser.c b/test/raytracer/gmlparser.c new file mode 100644 index 0000000..b0a6cbc --- /dev/null +++ b/test/raytracer/gmlparser.c @@ -0,0 +1,102 @@ +/* Parser for GML */ + +#include "config.h" + +#include "arrays.h" +#include "gml.h" +#include "gmllexer.h" +#include "gmlparser.h" + +static struct array * parse_tokenlist(void); + +static int parse_token(struct tok * t) +{ + get_lexeme(); + switch (current_lexeme.kind) { + case OPERATOR: + discard_lexeme(); + t->tag = current_lexeme.u.op; + break; + case IDENTIFIER: + discard_lexeme(); + t->tag = Identifier; + t->u.s = current_lexeme.u.s; + break; + case BINDER: + discard_lexeme(); + t->tag = Binder; + t->u.s = current_lexeme.u.s; + break; + case BOOLEAN: + discard_lexeme(); + t->tag = Boolean; + t->u.i = current_lexeme.u.i; + break; + case INTEGER: + discard_lexeme(); + t->tag = Integer; + t->u.i = current_lexeme.u.i; + break; + case REAL: + discard_lexeme(); + t->tag = Real; + t->u.d = current_lexeme.u.d; + break; + case STRING: + discard_lexeme(); + t->tag = String; + t->u.s = current_lexeme.u.s; + break; + case LBRACE: + discard_lexeme(); + t->tag = Function; + t->u.a = parse_tokenlist(); + get_lexeme(); + if (current_lexeme.kind != RBRACE) { + fprintf(stderr, "} expected\n"); + exit(2); + } + discard_lexeme(); + break; + case LBRACKET: + discard_lexeme(); + t->tag = Array; + t->u.a = parse_tokenlist(); + get_lexeme(); + if (current_lexeme.kind != RBRACKET) { + fprintf(stderr, "] expected\n"); + exit(2); + } + discard_lexeme(); + break; + default: + return 0; + } + return 1; +} + +static struct array * parse_tokenlist(void) +{ + struct array * a = alloc_array(sizeof(struct tok), 10); + struct tok t; + int i = 0; + while (parse_token(&t)) { + extend_array(struct tok, a); + set_array_large(struct tok, a, i, t); + i++; + } + return a; +} + +struct array * parse_program(void) +{ + struct array * a = parse_tokenlist(); + get_lexeme(); + if (current_lexeme.kind != END_OF_FILE) { + fprintf(stderr, "syntax error at end of program\n"); + exit(2); + } + return a; +} + + diff --git a/test/raytracer/gmlparser.h b/test/raytracer/gmlparser.h new file mode 100644 index 0000000..6a31546 --- /dev/null +++ b/test/raytracer/gmlparser.h @@ -0,0 +1,3 @@ +/* Parser for GML */ + +struct array * parse_program(void); diff --git a/test/raytracer/intersect.c b/test/raytracer/intersect.c new file mode 100644 index 0000000..c8a2126 --- /dev/null +++ b/test/raytracer/intersect.c @@ -0,0 +1,416 @@ +#include "config.h" +#include "arrays.h" +#include "point.h" +#include "vector.h" +#include "eval.h" +#include "object.h" +#include "matrix.h" +#include "intersect.h" + +/* Operations on interval lists */ + +#define POS_INFTY 1e300 + +struct intervlist { + flt beg; + struct object * obeg; + flt end; + struct object * oend; + struct intervlist * next; +}; + +typedef struct intervlist * intervlist; + +static intervlist cons(flt b, struct object * ob, + flt e, struct object * oe, + intervlist next) +{ + intervlist l = arena_alloc(sizeof(struct intervlist)); + l->beg = b; + l->obeg = ob; + l->end = e; + l->oend = oe; + l->next = next; + return l; +} + +static intervlist conshead(intervlist i, intervlist next) +{ + intervlist l = arena_alloc(sizeof(struct intervlist)); + l->beg = i->beg; + l->obeg = i->obeg; + l->end = i->end; + l->oend = i->oend; + l->next = next; + return l; +} + +static intervlist union_intervals(intervlist l1, intervlist l2) +{ + flt beg; + struct object * obeg; + + if (l1 == NULL) return l2; + if (l2 == NULL) return l1; + if (l1->end < l2->beg) + /* int1 strictly before int2 */ + return conshead(l1, union_intervals(l1->next, l2)); + if (l2->end < l1->beg) + /* int2 strictly before int1 */ + return conshead(l2, union_intervals(l1, l2->next)); + /* int1 and int2 overlap, merge */ + if (l1->beg < l2->beg) { + beg = l1->beg; + obeg = l1->obeg; + } else { + beg = l2->beg; + obeg = l2->obeg; + } + if (l1->end < l2->end) + return union_intervals(l1->next, + cons(beg, obeg, l2->end, l2->oend, l2->next)); + else + return union_intervals(cons(beg, obeg, l1->end, l1->oend, l1->next), + l2->next); +} + +static intervlist intersect_intervals(intervlist l1, intervlist l2) +{ + flt beg; + struct object * obeg; + + if (l1 == NULL) return NULL; + if (l2 == NULL) return NULL; + if (l1->end < l2->beg) + /* int1 strictly before int2 */ + return intersect_intervals(l1->next, l2); + if (l2->end < l1->beg) + /* int2 strictly before int1 */ + return intersect_intervals(l1, l2->next); + /* int1 and int2 overlap, add intersection */ + if (l1->beg > l2->beg) { + beg = l1->beg; + obeg = l1->obeg; + } else { + beg = l2->beg; + obeg = l2->obeg; + } + if (l1->end < l2->end) + return cons(beg, obeg, l1->end, l1->oend, intersect_intervals(l1->next, l2)); + else + return cons(beg, obeg, l2->end, l2->oend, intersect_intervals(l1, l2->next)); +} + +static intervlist difference_intervals(intervlist l1, intervlist l2) +{ + intervlist l; + + if (l1 == NULL) return NULL; + if (l2 == NULL) return l1; + if (l1->end < l2->beg) + /* int1 strictly before int2, keep int1 */ + return conshead(l1, difference_intervals(l1->next, l2)); + if (l2->end < l1->beg) + /* int2 strictly before int1, throw int2 away */ + return difference_intervals(l1, l2->next); + /* int and int2 overlap */ + if (l1->end > l2->end) + /* int1 extends beyond int2, add segment [l2->end,l1->end] */ + l = difference_intervals(cons(l2->end, l2->oend, l1->end, l1->oend, l1->next), l2->next); + else + /* int2 doesn't extend beyond int2, nothing to add */ + l = difference_intervals(l1->next, l2); + if (l1->beg < l2->beg) + /* int1 starts before l2, add segment [l1->beg,l2->beg] */ + l = cons(l1->beg, l1->obeg, l2->beg, l2->obeg, l); + return l; +} + +/* Intersections between rays (p + tv) and basic objects (bobj) */ + +static intervlist intersect_ray_plane(struct object * bobj, + struct point * p, + struct vector * v) +{ + if (p->y >= 0) + if (v->dy >= 0) + return NULL; + else + return cons(- p->y / v->dy, bobj, POS_INFTY, bobj, NULL); + else + if (v->dy >= 0) + return cons(0, bobj, - p->y / v->dy, bobj, NULL); + else + return cons(0, bobj, POS_INFTY, bobj, NULL); +} + +static intervlist intersect_ray_sphere(struct object * bobj, + struct point * p, + struct vector * v) +{ + flt a, b, c, d, sq, t0, t1, tswap; + + a = v->dx * v->dx + v->dy * v->dy + v->dz * v->dz; + b = p->x * v->dx + p->y * v->dy + p->z * v->dz; + c = p->x * p->x + p->y * p->y + p->z * p->z - 1.0; + d = b * b - a * c; + if (d <= 0) return NULL; + sq = sqrt(d); + /* Numerically stable solution of quadratic */ + if (b >= 0) { + t0 = c / (- b - sq); + t1 = (- b - sq) / a; + } else { + t0 = c / (- b + sq); + t1 = (- b + sq) / a; + } + if (t0 >= t1) { + tswap = t0; t0 = t1; t1 = tswap; + } + if (t1 <= 0) { + assert (t0 <= 0); + return NULL; + } + if (t0 >= 0) { + assert (t1 >= 0); + return cons(t0, bobj, t1, bobj, NULL); + } + return cons(0, bobj, t1, bobj, NULL); +} + +static intervlist intersect_ray_slice(struct object * bobj, + flt pc, flt vc) +{ + if (pc > 1.0) { + if (vc >= 0.0) + return NULL; + else + return cons((1.0 - pc) / vc, bobj, - pc / vc, bobj, NULL); + } + if (pc < 0.0) { + if (vc <= 0.0) + return NULL; + else + return cons(- pc / vc, bobj, (1.0 - pc) / vc, bobj, NULL); + } + if (vc == 0.0) + return cons(0.0, bobj, POS_INFTY, bobj, NULL); + if (vc > 0.0) + return cons(0.0, bobj, (1.0 - pc) / vc, bobj, NULL); + else + return cons(0.0, bobj, - pc / vc, bobj, NULL); +} + +static intervlist intersect_ray_cube(struct object * bobj, + struct point * p, + struct vector * v) +{ + return intersect_intervals(intersect_ray_slice(bobj, p->x, v->dx), + intersect_intervals(intersect_ray_slice(bobj, p->y, v->dy), + intersect_ray_slice(bobj, p->z, v->dz))); +} + +static intervlist intersect_ray_infinite_cylinder(struct object * bobj, + struct point * p, + struct vector * v) +{ + flt a, b, c, d, sq, t0, t1, tswap; + + a = v->dx * v->dx + v->dz * v->dz; + b = p->x * v->dx + p->z * v->dz; + c = p->x * p->x + p->z * p->z - 1.0; + d = b * b - a * c; + if (d <= 0.0) return NULL; + sq = sqrt(d); + if (b >= 0.0) { + t0 = c / (- b - sq); + t1 = (- b - sq) / a; + } else { + t0 = c / (- b + sq); + t1 = (- b + sq) / a; + } + if (t0 >= t1) { tswap = t0; t0 = t1; t1 = tswap; } + if (t1 <= 0.0) { + assert (t0 <= 0.0); + return NULL; + } + if (t0 >= 0.0) { + assert (t1 >= 0.0); + return cons(t0, bobj, t1, bobj, NULL); + } + return cons(0.0, bobj, t1, bobj, NULL); +} + +static intervlist intersect_ray_cylinder(struct object * bobj, + struct point * p, + struct vector * v) +{ + return intersect_intervals(intersect_ray_infinite_cylinder(bobj, p, v), + intersect_ray_slice(bobj, p->y, v->dy)); +} + +static intervlist intersect_ray_infinite_cone(struct object * bobj, + struct point * p, + struct vector * v) +{ + flt a, b, c, d, sq, t, t0, t1, tswap; + + a = v->dx * v->dx - v->dy * v->dy + v->dz * v->dz; + b = p->x * v->dx - p->y * v->dy + p->z * v->dz; + c = p->x * p->x - p->y * p->y + p->z * p->z; + if (a == 0.0) { + if (b == 0.0) return NULL; + t = - 0.5 * c / b; + if (c < 0.0) { + if (t <= 0.0) + return cons(0.0, bobj, POS_INFTY, bobj, NULL); + else + return cons(0.0, bobj, t, bobj, NULL); + } + if (t <= 0.0) + return NULL; + else + return cons(t, bobj, POS_INFTY, bobj, NULL); + } + d = b * b - a * c; + if (d <= 0.0) return NULL; + sq = sqrt(d); + if (b >= 0.0) { + t0 = c / (- b - sq); + t1 = (- b - sq) / a; + } else { + t0 = c / (- b + sq); + t1 = (- b + sq) / a; + } + if (t0 >= t1) { tswap = t0; t0 = t1; t1 = tswap; } + if (t1 <= 0.0) { + assert (t0 <= 0.0); + if (c < 0.0) + return cons(0.0, bobj, POS_INFTY, bobj, NULL); + else + return NULL; + } + if (t0 >= 0.0) { + assert (t1 >= 0.0); + if (c < 0.0) + return cons(0.0, bobj, t0, bobj, + cons (t1, bobj, POS_INFTY, bobj, NULL)); + else + return cons(t0, bobj, t1, bobj, NULL); + } + if (c < 0.0) + return cons(0.0, bobj, t1, bobj, NULL); + else + return cons(t1, bobj, POS_INFTY, bobj, NULL); +} + +static intervlist intersect_ray_cone(struct object * bobj, + struct point * p, + struct vector * v) +{ + return intersect_intervals(intersect_ray_infinite_cone(bobj, p, v), + intersect_ray_slice(bobj, p->y, v->dy)); +} + +/* Approximate test based on bounding sphere */ + +static inline int inside_bs(struct point * p, + struct vector * v, + struct object * obj) +{ + flt x, y, z, a, b, c; + + assert(obj->radius >= 0.0); + + if (obj->radius >= POS_INFTY) return 1; + /* Shift origin to obj.center */ + x = p->x - obj->center.x; + y = p->y - obj->center.y; + z = p->z - obj->center.z; + /* Check whether quadratic has positive discriminant */ + a = v->dx * v->dx + v->dy * v->dy + v->dz * v->dz; + b = x * v->dx + y * v->dy + z * v->dz; + c = x * x + y * y + z * z - obj->radius * obj->radius; + return (b * b > a * c); +} + +/* Interval list representing the intersection of the given ray + and the given composite object */ + +static intervlist intersect_ray_object(struct point * p, + struct vector * v, + struct object * obj) +{ +#define CONVERT_V_P \ + apply_to_point(obj->world2obj, p, &p2); \ + apply_to_vect(obj->world2obj, v, &v2) + struct point p2; + struct vector v2; + intervlist res; + + /* Fast, approximate test based on bounding sphere */ + if (! inside_bs(p, v, obj)) return NULL; + /* Slow, exact test */ + switch (obj->kind) { + case Cone: + CONVERT_V_P; + res = intersect_ray_cone(obj, &p2, &v2); + break; + case Cube: + CONVERT_V_P; + res = intersect_ray_cube(obj, &p2, &v2); + break; + case Cylinder: + CONVERT_V_P; + res = intersect_ray_cylinder(obj, &p2, &v2); + break; + case Plane: + CONVERT_V_P; + res = intersect_ray_plane(obj, &p2, &v2); + break; + case Sphere: + CONVERT_V_P; + res = intersect_ray_sphere(obj, &p2, &v2); + break; + case Union: + res = union_intervals(intersect_ray_object(p, v, obj->o1), + intersect_ray_object(p, v, obj->o2)); + break; + case Intersection: + res = intersect_intervals(intersect_ray_object(p, v, obj->o1), + intersect_ray_object(p, v, obj->o2)); + break; + case Difference: + res = difference_intervals(intersect_ray_object(p, v, obj->o1), + intersect_ray_object(p, v, obj->o2)); + break; + default: + assert(0); + } +#undef CONVERT_V_P + return res; +} + +/* Return the closest base object intersecting the given ray, and + the curvilinear abscissa t of the intersection point. + Return NULL if no intersection. + Return t = 0.0 if the viewpoint (origin of the ray) is inside an + object. */ + +struct object * intersect_ray(struct point * p, + struct vector * v, + struct object * obj, + int initial, + /*out*/ flt * t) +{ + intervlist l = intersect_ray_object(p, v, obj); + if (l == NULL) return NULL; + assert(l->beg >= 0.0); + if (l->beg <= 0.0 && !initial) { + /* skip first intersection */ + l = l->next; + if (l == NULL) return NULL; + } + *t = l->beg; + return l->obeg; +} diff --git a/test/raytracer/intersect.h b/test/raytracer/intersect.h new file mode 100644 index 0000000..35241a6 --- /dev/null +++ b/test/raytracer/intersect.h @@ -0,0 +1,11 @@ +/* Return the closest base object intersecting the given ray, and + the curvilinear abscissa t of the intersection point. + Return NULL if no intersection. + Return t = 0.0 if the viewpoint (origin of the ray) is inside an + object. */ + +struct object * intersect_ray(struct point * p, + struct vector * v, + struct object * obj, + int initial, + /*out*/ flt * t); diff --git a/test/raytracer/kal.gml b/test/raytracer/kal.gml new file mode 100644 index 0000000..e9d5815 --- /dev/null +++ b/test/raytracer/kal.gml @@ -0,0 +1,78 @@ +% kal.gml + +1.0 0.7 0.7 point /red +0.7 1.0 0.7 point /green +0.7 0.7 1.0 point /blue + +{ /color + { /v /u /face + color 0.1 0.99 1.0 + } +} /mirror + + % a cube +{ /v /u /face % bind arguments + 0.9 1.0 0.9 point % surface color + 0.1 0.9 1.0 % kd ks n +} cube + +0.0 -0.5 35.0 translate + % a sphere +{ /v /u /face % bind arguments + 0.9 0.9 1.0 point % surface color + 0.1 0.9 1.0 % kd ks n +} sphere + +-1.0 0.0 20.0 translate + +union + % a sphere +{ /v /u /face % bind arguments + 1.0 1.0 0.9 point % surface color + 0.1 0.9 1.0 % kd ks n +} cylinder +20.0 rotatex +0.5 0.5 15.0 translate +union + +blue mirror apply plane +0.0 -2.0 0.0 translate + +union + +green mirror apply plane +0.0 -2.0 0.0 translate +120.0 rotatez + +union + +red mirror apply plane +0.0 -2.0 0.0 translate +240.0 rotatez + +union + +{ /v /u /face + 0.4 0.4 0.4 point + 0.0 0.0 1.0 +} plane + +90.0 rotatex +0.0 0.0 3.0 translate +difference + +0.0 -1.0 0.0 translate + +/scene + +0.0 0.0 -2.0 point 0.9 0.9 0.9 point pointlight /light1 + +0.4 0.4 0.4 point % ambient +[light1] % lights +scene % object +1000 % depth +90.0 % fov +400 400 % wid ht +"kal.ppm" % output file +render + diff --git a/test/raytracer/light.c b/test/raytracer/light.c new file mode 100644 index 0000000..1c90104 --- /dev/null +++ b/test/raytracer/light.c @@ -0,0 +1,126 @@ +#include "config.h" +#include "point.h" +#include "vector.h" +#include "eval.h" +#include "object.h" +#include "intersect.h" +#include "light.h" + +struct light * dirlight(struct point * dir, struct point * c) +{ + struct light * l = arena_alloc(sizeof(struct light)); + l->kind = Directional; + ASSIGN(l->u.directional.dir, *dir); + ASSIGN(l->u.directional.col, *c); + return l; +} + +struct light * pointlight(struct point * orig, struct point * c) +{ + struct light * l = arena_alloc(sizeof(struct light)); + l->kind = Pointlight; + ASSIGN(l->u.point.orig, *orig); + ASSIGN(l->u.point.col, *c); + return l; +} + +struct light * spotlight(struct point * pos, struct point * at, + struct point * c, flt cutoff, flt exp) +{ + struct light * l = arena_alloc(sizeof(struct light)); + struct vector uv; + l->kind = Spot; + ASSIGN(l->u.spot.orig, *pos); + ASSIGN(l->u.spot.at, *at); + ASSIGN(l->u.spot.col, *c); + l->u.spot.cutoff = cutoff; + l->u.spot.exponent = exp; + between(at, pos, &uv); + vnormalize(&uv, &(l->u.spot.unit)); + return l; +} + +static void color_from_light(struct object * scene, + struct point * pt, + struct vector * v, + struct vector * n, + flt kd, flt ks, flt phong, + struct light * light, + /*out*/ struct point * il) +{ + struct vector dir; + struct object * obj; + struct vector L, H, vnorm; + struct point i; + flt t, att, dp, m; + + /* Compute direction towards light source */ + switch (light->kind) { + case Directional: + dir.dx = - light->u.directional.dir.x; + dir.dy = - light->u.directional.dir.y; + dir.dz = - light->u.directional.dir.z; + break; + case Pointlight: + between(pt, &light->u.point.orig, &dir); + break; + case Spot: + between(pt, &light->u.spot.orig, &dir); + break; + } + /* Check that light source is not behind us */ + if (dotproduct(&dir, n) <= 0.0) return; /*no contribution*/ + /* Check that no object blocks the light */ + obj = intersect_ray(pt, &dir, scene, 0 /*??*/, &t); + if (obj != NULL && (light->kind == Directional || t <= 1.0)) return; + /* Compute the L and H vectors */ + vnormalize(&dir, &L); + vnormalize(v, &vnorm); + vsub(&L, &vnorm, &H); + vnormalize(&H, &H); + /* Intensity of light source at object */ + switch (light->kind) { + case Directional: + i.x = light->u.directional.col.x; + i.y = light->u.directional.col.y; + i.z = light->u.directional.col.z; + break; + case Pointlight: + att = 100.0 / (99.0 + dist2(pt, &light->u.point.orig)); + i.x = light->u.point.col.x * att; + i.y = light->u.point.col.y * att; + i.z = light->u.point.col.z * att; + break; + case Spot: + dp = dotproduct(&L, &light->u.spot.unit); + if (acos(dp) >= light->u.spot.cutoff) return; + att = + pow(dp, light->u.spot.exponent) + * (100.0 / (99.0 + dist2(pt, &light->u.spot.orig))); + i.x = light->u.spot.col.x * att; + i.y = light->u.spot.col.y * att; + i.z = light->u.spot.col.z * att; + break; + } + /* Add final contribution */ + m = kd * dotproduct(n, &L) + ks * pow(dotproduct(n, &H), phong); + il->x += m * i.x; + il->y += m * i.y; + il->z += m * i.z; +} + +void color_from_lights(struct object * scene, + struct point * p, + struct vector * v, + struct vector * n, + flt kd, flt ks, flt phong, + struct light ** lights, + /*out*/ struct point * il) +{ + int i; + + il->x = il->y = il->z = 0.0; + for (i = 0; lights[i] != NULL; i++) + color_from_light(scene, p, v, n, kd, ks, phong, lights[i], il); +} + diff --git a/test/raytracer/light.h b/test/raytracer/light.h new file mode 100644 index 0000000..9a1d3e2 --- /dev/null +++ b/test/raytracer/light.h @@ -0,0 +1,33 @@ +struct light { + enum { Directional, Pointlight, Spot } kind; + union { + struct { + struct point dir; + struct point col; + } directional; + struct { + struct point orig; + struct point col; + } point; + struct { + struct point orig; + struct point at; + struct point col; + flt cutoff, exponent; + struct vector unit; + } spot; + } u; +}; + +struct light * dirlight(struct point * dir, struct point * c); +struct light * pointlight(struct point * orig, struct point * c); +struct light * spotlight(struct point * pos, struct point * at, + struct point * c, flt a, flt cutoff); + +void color_from_lights(struct object * obj, + struct point * p, + struct vector * v, + struct vector * n, + flt kd, flt ks, flt phong, + struct light ** lights, + /*out*/ struct point * il); diff --git a/test/raytracer/main.c b/test/raytracer/main.c new file mode 100644 index 0000000..1efd396 --- /dev/null +++ b/test/raytracer/main.c @@ -0,0 +1,13 @@ +#include "config.h" +#include "arrays.h" +#include "gmllexer.h" +#include "gmlparser.h" +#include "eval.h" + +int main(int argc, char ** argv) +{ + arena_init(); + init_lexer(); + execute_program(parse_program()); + return 0; +} diff --git a/test/raytracer/matrix.c b/test/raytracer/matrix.c new file mode 100644 index 0000000..8814346 --- /dev/null +++ b/test/raytracer/matrix.c @@ -0,0 +1,102 @@ +#include "config.h" +#include "point.h" +#include "vector.h" +#include "matrix.h" + +struct matrix matrix_identity = { + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, +}; + +void apply_to_point(struct matrix * m, struct point * p, + /*out*/ struct point * r) +{ + r->x = m->xx * p->x + m->xy * p->y + m->xz * p->z + m->xt; + r->y = m->yx * p->x + m->yy * p->y + m->yz * p->z + m->yt; + r->z = m->zx * p->x + m->zy * p->y + m->zz * p->z + m->zt; +} + +void apply_to_vect(struct matrix * m, struct vector * v, + /*out*/ struct vector * r) +{ + r->dx = m->xx * v->dx + m->xy * v->dy + m->xz * v->dz; + r->dy = m->yx * v->dx + m->yy * v->dy + m->yz * v->dz; + r->dz = m->zx * v->dx + m->zy * v->dy + m->zz * v->dz; +} + + + +struct matrix * mtranslate(flt sx, flt sy, flt sz) +{ + struct matrix * m = arena_alloc(sizeof(struct matrix)); + m->xx = 1.0; m->xy = 0.0; m->xz = 0.0; m->xt = sx; + m->yx = 0.0; m->yy = 1.0; m->yz = 0.0; m->yt = sy; + m->zx = 0.0; m->zy = 0.0; m->zz = 1.0; m->zt = sz; + return m; +} + +struct matrix * mscale(flt sx, flt sy, flt sz) +{ + struct matrix * m = arena_alloc(sizeof(struct matrix)); + m->xx = sx; m->xy = 0.0; m->xz = 0.0; m->xt = 0.0; + m->yx = 0.0; m->yy = sy ; m->yz = 0.0; m->yt = 0.0; + m->zx = 0.0; m->zy = 0.0; m->zz = sz ; m->zt = 0.0; + return m; +} + +struct matrix * mrotatex(flt a) +{ + struct matrix * m = arena_alloc(sizeof(struct matrix)); + flt c = cos(a); + flt s = sin(a); + m->xx = 1.0; m->xy = 0.0; m->xz = 0.0; m->xt = 0.0; + m->yx = 0.0; m->yy = c; m->yz = - s; m->yt = 0.0; + m->zx = 0.0; m->zy = s; m->zz = c; m->zt = 0.0; + return m; +} + +struct matrix * mrotatey(flt a) +{ + struct matrix * m = arena_alloc(sizeof(struct matrix)); + flt c = cos(a); + flt s = sin(a); + m->xx = c; m->xy = 0.0; m->xz = s; m->xt = 0.0; + m->yx = 0.0; m->yy = 1.0; m->yz = 0.0; m->yt = 0.0; + m->zx = - s; m->zy = 0.0; m->zz = c; m->zt = 0.0; + return m; +} + +struct matrix * mrotatez(flt a) +{ + struct matrix * m = arena_alloc(sizeof(struct matrix)); + flt c = cos(a); + flt s = sin(a); + m->xx = c; m->xy = - s; m->xz = 0.0; m->xt = 0.0; + m->yx = s; m->yy = c; m->yz = 0.0; m->yt = 0.0; + m->zx = 0.0; m->zy = 0.0; m->zz = 1.0; m->zt = 0.0; + return m; +} + +struct matrix * mcompose(struct matrix * m, struct matrix * n) +{ + struct matrix * r = arena_alloc(sizeof(struct matrix)); + + r->xx = m->xx * n->xx + m->xy * n->yx + m->xz * n->zx; + r->xy = m->xx * n->xy + m->xy * n->yy + m->xz * n->zy; + r->xz = m->xx * n->xz + m->xy * n->yz + m->xz * n->zz; + r->xt = m->xx * n->xt + m->xy * n->yt + m->xz * n->zt + m->xt; + + r->yx = m->yx * n->xx + m->yy * n->yx + m->yz * n->zx; + r->yy = m->yx * n->xy + m->yy * n->yy + m->yz * n->zy; + r->yz = m->yx * n->xz + m->yy * n->yz + m->yz * n->zz; + r->yt = m->yx * n->xt + m->yy * n->yt + m->yz * n->zt + m->yt; + + r->zx = m->zx * n->xx + m->zy * n->yx + m->zz * n->zx; + r->zy = m->zx * n->xy + m->zy * n->yy + m->zz * n->zy; + r->zz = m->zx * n->xz + m->zy * n->yz + m->zz * n->zz; + r->zt = m->zx * n->xt + m->zy * n->yt + m->zz * n->zt + m->zt; + + return r; +} + diff --git a/test/raytracer/matrix.h b/test/raytracer/matrix.h new file mode 100644 index 0000000..4b9f434 --- /dev/null +++ b/test/raytracer/matrix.h @@ -0,0 +1,18 @@ +struct matrix { + flt xx, xy, xz, xt; + flt yx, yy, yz, yt; + flt zx, zy, zz, zt; +}; + +void apply_to_point(struct matrix * m, struct point * p, + /*out*/ struct point * r); +void apply_to_vect(struct matrix * m, struct vector * v, + /*out*/ struct vector * r); +extern struct matrix matrix_identity; +#define mid (&matrix_identity) +struct matrix * mtranslate(flt sx, flt sy, flt sz); +struct matrix * mscale(flt sx, flt sy, flt sz); +struct matrix * mrotatex(flt a); +struct matrix * mrotatey(flt a); +struct matrix * mrotatez(flt a); +struct matrix * mcompose(struct matrix * m, struct matrix * n); diff --git a/test/raytracer/memory.c b/test/raytracer/memory.c new file mode 100644 index 0000000..fe194cc --- /dev/null +++ b/test/raytracer/memory.c @@ -0,0 +1,60 @@ +/* Arena-based memory management */ + +#include "config.h" + +#define SIZE_AREA 1024000 + +struct area { + struct area * next; + char data[SIZE_AREA]; +}; + +struct area * arena_head; +struct area * arena_cur; +int arena_cur_ofs; + +void arena_init(void) +{ + arena_head = malloc(sizeof(struct area)); + if (arena_head == NULL) { + fprintf(stderr, "Out of memory\n"); + exit(2); + } + arena_head->next = NULL; + arena_cur = arena_head; + arena_cur_ofs = 0; +} + +void arena_clear(void) +{ + arena_cur = arena_head; + arena_cur_ofs = 0; +} + +void * arena_alloc(int size) +{ + char * res; + struct area * n; + + if (size >= SIZE_AREA) { + fprintf(stderr, "Memory request too large (%d)\n", size); + exit(2); + } + if (arena_cur_ofs + size <= SIZE_AREA) { + res = arena_cur->data + arena_cur_ofs; + arena_cur_ofs += size; + return res; + } + if (arena_cur->next == NULL) { + n = malloc(sizeof(struct area)); + if (n == NULL) { + fprintf(stderr, "Out of memory\n"); + exit(2); + } + n->next = NULL; + arena_cur->next = n; + } + arena_cur = arena_cur->next; + arena_cur_ofs = size; + return arena_cur->data; +} diff --git a/test/raytracer/object.c b/test/raytracer/object.c new file mode 100644 index 0000000..3dd6e77 --- /dev/null +++ b/test/raytracer/object.c @@ -0,0 +1,214 @@ +#include "config.h" +#include "arrays.h" +#include "point.h" +#include "vector.h" +#include "eval.h" +#include "object.h" +#include "matrix.h" + +static struct object * new_object(int kind, struct closure * c) +{ + struct object * o = arena_alloc(sizeof(struct object)); + o->kind = kind; + ASSIGN(o->surf, *c); + o->world2obj = o->obj2world = mid; + o->max_scale_applied = 1.0; + o->radius = BS_NOT_COMPUTED; + return o; +} + +struct object * cone(struct closure * c) +{ return new_object(Cone, c); } + +struct object * cube(struct closure * c) +{ return new_object(Cube, c); } + +struct object * cylinder(struct closure * c) +{ return new_object(Cylinder, c); } + +struct object * plane(struct closure * c) +{ return new_object(Plane, c); } + +struct object * sphere(struct closure * c) +{ return new_object(Sphere, c); } + +static struct object * transform(struct object * o, + struct matrix * t, + struct matrix * tinv, + flt scale) +{ + struct object * no = arena_alloc(sizeof(struct object)); + no->kind = o->kind; + switch (o->kind) { + case Union: + case Intersection: + case Difference: + no->o1 = transform(o->o1, t, tinv, scale); + no->o2 = transform(o->o2, t, tinv, scale); + break; + default: + ASSIGN(no->surf, o->surf); + no->world2obj = mcompose(o->world2obj, tinv); + no->obj2world = mcompose(t, o->obj2world); + no->max_scale_applied = o->max_scale_applied * scale; + } + no->radius = BS_NOT_COMPUTED; + return no; +} + +struct object * orotatex(struct object * o1, flt a) +{ + return transform(o1, mrotatex(a), mrotatex(-a), 1.0); +} + +struct object * orotatey(struct object * o1, flt a) +{ + return transform(o1, mrotatey(a), mrotatey(-a), 1.0); +} + +struct object * orotatez(struct object * o1, flt a) +{ + return transform(o1, mrotatez(a), mrotatez(-a), 1.0); +} + +struct object * oscale(struct object * o1, flt sx, flt sy, flt sz) +{ + flt scale = sx; + if (sy > scale) scale = sy; + if (sz > scale) scale = sz; + return transform(o1, mscale(sx, sy, sz), mscale(1 / sx, 1 / sy, 1 / sz), + scale); +} + +struct object * otranslate(struct object * o1, + flt tx, flt ty, flt tz) +{ + return transform(o1, mtranslate(tx, ty, tz), mtranslate(- tx, - ty, - tz), + 1.0); +} + +struct object * ouscale(struct object * o1, flt s) +{ + flt sinv = 1 / s; + return transform(o1, mscale(s, s, s), mscale(sinv, sinv, sinv), s); +} + +struct object * odifference(struct object * o1, struct object * o2) +{ + struct object * o = arena_alloc(sizeof(struct object)); + o->kind = Difference; + o->o1 = o1; + o->o2 = o2; + o->radius = BS_NOT_COMPUTED; + return o; +} + +struct object * ointersect(struct object * o1, struct object * o2) +{ + struct object * o = arena_alloc(sizeof(struct object)); + o->kind = Intersection; + o->o1 = o1; + o->o2 = o2; + o->radius = BS_NOT_COMPUTED; + return o; +} + +struct object * ounion(struct object * o1, struct object * o2) +{ + struct object * o = arena_alloc(sizeof(struct object)); + o->kind = Union; + o->o1 = o1; + o->o2 = o2; + o->radius = BS_NOT_COMPUTED; + return o; +} + +static void normal_vector_object(struct object * obj, + struct point * p, + int face, + /*out*/ struct vector * v) +{ + switch (obj->kind) { + case Cone: + if (face == 0) { + v->dx = p->x; v->dy = - p->y; v->dz = p->z; + } else { + v->dx = 0; v->dy = 1; v->dz = 0; + } + break; + case Cube: + switch (face) { + case 0: + v->dx = 0; v->dy = 0; v->dz = -1; break; + case 1: + v->dx = 0; v->dy = 0; v->dz = 1; break; + case 2: + v->dx = -1; v->dy = 0; v->dz = 0; break; + case 3: + v->dx = 1; v->dy = 0; v->dz = 0; break; + case 4: + v->dx = 0; v->dy = 1; v->dz = 0; break; + case 5: + v->dx = 0; v->dy = -1; v->dz = 0; break; + } + break; + case Cylinder: + switch (face) { + case 0: + v->dx = p->x; v->dy = 0; v->dz = p->z; break; + case 1: + v->dx = 0; v->dy = 1; v->dz = 0; break; + case 2: + v->dx = 0; v->dy = -1; v->dz = 0; break; + } + break; + case Plane: + v->dx = 0; v->dy = 1; v->dz = 0; break; + case Sphere: + v->dx = p->x; v->dy = p->y; v->dz = p->z; break; + default: + assert(0); + } +} + +static void tangent_vectors(struct vector * n, + /*out*/ struct vector * t1, + /*out*/ struct vector * t2) +{ + if (n->dy > 0) { + t1->dx = n->dy; t1->dy = - n->dx; t1->dz = 0; + t2->dx = 0; t2->dy = n->dz; t2->dz = - n->dy; + } + /* y 0 xy x + -x ^ z = yy = y y + 0 -y yz z */ + else if (n->dy == 0) { + t1->dx = n->dz; t1->dy = 0; t1->dz = - n->dx; + t2->dx = n->dz; t2->dy = 1; t2->dz = - n->dx; + } + /* z z x x + 0 ^ 1 = 0 = 1 y + -x -x z z */ + else { + t1->dx = n->dy; t1->dy = - n->dx; t1->dz = 0; + t2->dx = 0; t2->dy = - n->dz; t2->dz = n->dy; + } + /* y 0 -xy x + -x ^ -z = -yy = (-y) y + 0 y -zy z */ +} + +void normal_vector(struct object * obj, struct point * p, int face, + /*out*/ struct vector * n) +{ + struct point pobj; + struct vector nobj, tang_obj1, tang_obj2, tang_world1, tang_world2; + + apply_to_point(obj->world2obj, p, &pobj); + normal_vector_object(obj, &pobj, face, &nobj); + tangent_vectors(&nobj, &tang_obj1, &tang_obj2); + apply_to_vect(obj->obj2world, &tang_obj1, &tang_world1); + apply_to_vect(obj->obj2world, &tang_obj2, &tang_world2); + product(&tang_world1, &tang_world2, n); + vnormalize(n, n); +} diff --git a/test/raytracer/object.h b/test/raytracer/object.h new file mode 100644 index 0000000..d890bda --- /dev/null +++ b/test/raytracer/object.h @@ -0,0 +1,36 @@ +struct object { + enum { Cone, Cube, Cylinder, Plane, Sphere, + Union, Intersection, Difference } kind; + /* For base objects: Cone, Cube, Cylinder, Plane, Sphere */ + struct closure surf; /* surface function */ + struct matrix * world2obj, * obj2world; + flt max_scale_applied; + /* For compound objects: Union, Intersection, Difference */ + struct object * o1, * o2; + /* For all objects: bounding sphere (center and radius) */ + struct point center; + flt radius; +}; + +struct object * cone(struct closure * c); +struct object * cube(struct closure * c); +struct object * cylinder(struct closure * c); +struct object * plane(struct closure * c); +struct object * sphere(struct closure * c); + +struct object * orotatex(struct object * o1, flt a); +struct object * orotatey(struct object * o1, flt a); +struct object * orotatez(struct object * o1, flt a); +struct object * oscale(struct object * o1, flt sx, flt sy, flt sz); +struct object * otranslate(struct object * o1, + flt tx, flt ty, flt tz); +struct object * ouscale(struct object * o1, flt s); + +struct object * odifference(struct object * o1, struct object * o2); +struct object * ointersect(struct object * o1, struct object * o2); +struct object * ounion(struct object * o1, struct object * o2); + +void normal_vector(struct object * obj, struct point * p, int face, + /*out*/ struct vector * n); + +#define BS_NOT_COMPUTED (-1.0) diff --git a/test/raytracer/point.h b/test/raytracer/point.h new file mode 100644 index 0000000..b8d6ff0 --- /dev/null +++ b/test/raytracer/point.h @@ -0,0 +1,18 @@ +struct point { + flt x, y, z; +}; + +#if 0 +static inline flt dist2(struct point * p1, struct point * p2) +{ + flt dx = p2->x - p1->x; + flt dy = p2->y - p1->y; + flt dz = p2->z - p1->z; + return dx * dx + dy * dy + dz * dz; +} +#else +#define dist2(p1,p2) \ + (((p2)->x - (p1)->x) * ((p2)->x - (p1)->x) + \ + ((p2)->y - (p1)->y) * ((p2)->y - (p1)->y) + \ + ((p2)->z - (p1)->z) * ((p2)->z - (p1)->z)) +#endif diff --git a/test/raytracer/render.c b/test/raytracer/render.c new file mode 100644 index 0000000..3678cb6 --- /dev/null +++ b/test/raytracer/render.c @@ -0,0 +1,128 @@ +#include "config.h" +#include "point.h" +#include "vector.h" +#include "matrix.h" +#include "eval.h" +#include "object.h" +#include "light.h" +#include "intersect.h" +#include "surface.h" +#include "simplify.h" +#include "render.h" + +static void render_ray(struct point * amb, + int depth, + int initial, + struct object * obj, + struct light ** lights, + struct point * p, + struct vector * v, + flt multx, + flt multy, + flt multz, + /*out*/ struct point * color) +{ + struct object * bobj; + flt t; + struct point inter_w, inter_o; + int face; + flt surf_u, surf_v; + struct surface_characteristics sc; + flt dotprod; + struct vector n, n2, s; + struct point il, is; + + if (depth < 0) { + color->x = color->y = color->z = 0.0; + return; + } + bobj = intersect_ray(p, v, obj, initial, &t); + if (bobj == NULL || t == 0.0) { + color->x = color->y = color->z = 0.0; + return; + } + /* Compute surface characteristics */ + point_along(p, v, t, &inter_w); + apply_to_point(bobj->world2obj, &inter_w, &inter_o); + surface_coords(bobj, &inter_o, &face, &surf_u, &surf_v); + surface_function(&bobj->surf, face, surf_u, surf_v, &sc); + /* Construct the vectors on figure 4 */ + normal_vector(bobj, &inter_w, face, &n); + dotprod = dotproduct(v, &n); + vscale(&n, 2.0 * dotprod, &n2); + vsub(v, &n2, &s); + if (dotprod > 0.0) opposite(&n, &n); + /* Light sources */ + color_from_lights(obj, &inter_w, v, &n, + sc.kd, sc.ks, sc.phong, + lights, &il); + /* Recursive call for ray along s */ + multx = multx * sc.ks * sc.x; + multy = multy * sc.ks * sc.y; + multz = multz * sc.ks * sc.z; + if (multx < 0.1 && multy < 0.1 && multz < 0.1) + is.x = is.y = is.z = 0.0; + else + render_ray(amb, depth - 1, 0, obj, lights, &inter_w, &s, + multx, multy, multz, &is); + /* Compute final color */ + color->x = (sc.kd * amb->x + il.x + sc.ks * is.x) * sc.x; + color->y = (sc.kd * amb->y + il.y + sc.ks * is.y) * sc.y; + color->z = (sc.kd * amb->z + il.z + sc.ks * is.z) * sc.z; +} + +static int convert_color(flt c) +{ + int n = (int) (c * 255.0); + if (n < 0) n = 0; + if (n > 255) n = 255; + return n; +} + +void render(struct point * amb, + int numlights, + struct light ** lights, + struct object * scene, + int depth, + flt fov, + int wid, + int ht, + char * filename) +{ + FILE * oc; + flt wid2, ht2, scale, x, y; + int i, j; + struct point p; + struct vector v; + struct point color; + char * command; + + compute_bounding_spheres(scene); + wid2 = (flt) wid / 2.0; + ht2 = (flt) ht / 2.0; + scale = tan(fov / 2.0) / wid2; + oc = fopen(filename, "w"); + fprintf(oc, "P6\n# Camls 'R Us\n%d %d\n255\n", wid, ht); + arena_init(); + for (j = ht - 1; j >= 0; j--) { + y = ((flt) j - ht2 + 0.5) * scale; + for (i = 0; i < wid; i++) { + x = ((flt) i - wid2 + 0.5) * scale; + p.x = p.y = 0.0; p.z = -1.0; + v.dx = x; v.dy = y; v.dz = 1.0; + render_ray(amb, depth, 1, scene, lights, &p, &v, 255.0, 255.0, 255.0, + &color); + fputc(convert_color(color.x), oc); + fputc(convert_color(color.y), oc); + fputc(convert_color(color.z), oc); + arena_clear(); + } + } + fclose(oc); +#ifdef XV + command = malloc(strlen(filename) + 20); + sprintf(command, "xv %s &", filename); + system(command); + free(command); +#endif +} diff --git a/test/raytracer/render.h b/test/raytracer/render.h new file mode 100644 index 0000000..558e252 --- /dev/null +++ b/test/raytracer/render.h @@ -0,0 +1,9 @@ +void render(struct point * amb, + int numlights, + struct light ** lights, + struct object * scene, + int depth, + flt fov, + int wid, + int ht, + char * filename); diff --git a/test/raytracer/simplify.c b/test/raytracer/simplify.c new file mode 100644 index 0000000..7a4a545 --- /dev/null +++ b/test/raytracer/simplify.c @@ -0,0 +1,177 @@ +#include "config.h" +#include "point.h" +#include "vector.h" +#include "arrays.h" +#include "matrix.h" +#include "eval.h" +#include "object.h" +#include "simplify.h" + +#define INFINITE_RADIUS 1e300 + +static flt cone_radius = 1.0; +static flt cube_radius = 0.86602540378443859659; /* sqrt(3)/2 */ +static flt cylinder_radius = 1.11803398874989490253; /* sqrt(5)/2 */ +static flt sphere_radius = 1.0; + +static struct point cone_center = { 0.0, 1.0, 0.0 }; +static struct point cube_center = { 0.5, 0.5, 0.5 }; +static struct point cylinder_center = { 0.0, 0.5, 0.0 }; +static struct point sphere_center = { 0, 0, 0 }; +static struct point plane_center = { 0, 0, 0 }; + +static struct point origin = { 0, 0, 0 }; + +static inline void set_infinite(struct object * t) +{ + t->radius = INFINITE_RADIUS; +} + +static inline void union_bs(struct object * t1, struct object * t2, + struct object * obj) +{ + struct vector c1c2; + flt dd2, rr, rr2, d, alpha; + + if (t1->radius >= INFINITE_RADIUS || t2->radius >= INFINITE_RADIUS) { + obj->radius = INFINITE_RADIUS; + return; + } + between(&t1->center, &t2->center, &c1c2); + dd2 = vlength2(&c1c2); + rr = t2->radius - t1->radius; + rr2 = rr * rr; + if (dd2 <= rr2) { + /* take the biggest sphere */ + if (t1->radius <= t2->radius) { + ASSIGN(obj->center, t2->center); + obj->radius = t2->radius; + set_infinite(t2); + } else { + ASSIGN(obj->center, t1->center); + obj->radius = t1->radius; + set_infinite(t1); + } + return; + } + d = sqrt(dd2); + alpha = rr / (2 * d) + 0.5; + point_along(&t1->center, &c1c2, alpha, &obj->center); + obj->radius = (d + t1->radius + t2->radius) / 2; +} + +static inline void intersection_bs(struct object * t1, struct object * t2, + struct object * obj) +{ + struct vector c1c2; + flt dd2, rr, rr2, rpr, rpr2, diff, d, te1, te2, te3, te4, te, alpha; + + if (t1->radius >= INFINITE_RADIUS) { + ASSIGN(obj->center, t2->center); + obj->radius = t2->radius; + return; + } + if (t2->radius >= INFINITE_RADIUS) { + ASSIGN(obj->center, t1->center); + obj->radius = t1->radius; + return; + } + between(&t1->center, &t2->center, &c1c2); + dd2 = vlength2(&c1c2); + rr = t1->radius - t2->radius; + rr2 = rr * rr; + if (dd2 <= rr2) { + /* take the smallest sphere */ + if (t2->radius <= t1->radius) { + ASSIGN(obj->center, t2->center); + obj->radius = t2->radius; + set_infinite(t2); + } else { + ASSIGN(obj->center, t1->center); + obj->radius = t1->radius; + set_infinite(t1); + } + return; + } + rpr = t1->radius + t2->radius; + rpr2 = rpr * rpr; + if (dd2 > rpr2) { + ASSIGN(obj->center, origin); + obj->radius = 0.0; + return; + } + diff = t1->radius * t1->radius - t2->radius * t2->radius; + if (dd2 <= diff) { + ASSIGN(obj->center, t2->center); + obj->radius = t2->radius; + set_infinite(t2); + return; + } + if (dd2 <= -diff) { + ASSIGN(obj->center, t1->center); + obj->radius = t1->radius; + set_infinite(t1); + return; + } + d = sqrt(dd2); + te1 = t1->radius + d + t2->radius; + te2 = t1->radius + d - t2->radius; + te3 = t2->radius + d - t1->radius; + te4 = t1->radius + t2->radius - d; + te = (sqrt (te1 * te2 * te3 * te4)) / (2 * d); + alpha = + (t1->radius * t1->radius - t2->radius * t2->radius) / (2 * dd2) + 0.5; + point_along(&t1->center, &c1c2, alpha, &obj->center); + obj->radius = te; +} + +static inline void difference_bs(struct object * t1, struct object * t2, + struct object * obj) +{ + ASSIGN(obj->center, t1->center); + obj->radius = t1->radius; + set_infinite(t1); +} + +void compute_bounding_spheres(struct object * obj) +{ + if (obj->radius >= 0.0) return; /* already computed */ + switch (obj->kind) { + case Cone: + apply_to_point(obj->obj2world, &cone_center, &obj->center); + obj->radius = obj->max_scale_applied * cone_radius; + break; + case Cube: + apply_to_point(obj->obj2world, &cube_center, &obj->center); + obj->radius = obj->max_scale_applied * cube_radius; + break; + case Cylinder: + apply_to_point(obj->obj2world, &cylinder_center, &obj->center); + obj->radius = obj->max_scale_applied * cylinder_radius; + break; + case Plane: + ASSIGN(obj->center, plane_center); + obj->radius = INFINITE_RADIUS; + break; + case Sphere: + apply_to_point(obj->obj2world, &sphere_center, &obj->center); + obj->radius = obj->max_scale_applied * sphere_radius; + break; + case Union: + compute_bounding_spheres(obj->o1); + compute_bounding_spheres(obj->o2); + union_bs(obj->o1, obj->o2, obj); + break; + case Intersection: + compute_bounding_spheres(obj->o1); + compute_bounding_spheres(obj->o2); + intersection_bs(obj->o1, obj->o2, obj); + break; + case Difference: + compute_bounding_spheres(obj->o1); + compute_bounding_spheres(obj->o2); + difference_bs(obj->o1, obj->o2, obj); + break; + } +} + diff --git a/test/raytracer/simplify.h b/test/raytracer/simplify.h new file mode 100644 index 0000000..b2d1795 --- /dev/null +++ b/test/raytracer/simplify.h @@ -0,0 +1,3 @@ +/* Add bounding sphere info to the given object (recursively) */ + +void compute_bounding_spheres(struct object * obj); diff --git a/test/raytracer/surface.c b/test/raytracer/surface.c new file mode 100644 index 0000000..0059db1 --- /dev/null +++ b/test/raytracer/surface.c @@ -0,0 +1,148 @@ +#include "config.h" +#include "point.h" +#include "vector.h" +#include "eval.h" +#include "object.h" +#include "surface.h" + +/* Sphere: + x = sqrt(1 - y^2) sin(360u) + y = 2v - 1 + z = sqrt(1 - y^2) cos(360u) + hence v = (y+1)/2 + and u = atan2_turns(x, z) (atan2 in "number of turns") +*/ + +#define INV2PI (1.0 / (2.0 * M_PI)) + +static inline flt atan2_turns(flt x, flt z) +{ + flt r = atan2(x, z) * INV2PI; + if (r < 0.0) r += 1.0; + return r; +} + +static inline void sphere_coords(flt x, flt y, flt z, + int * face, flt * u, flt * v) +{ + *face = 0; + *u = atan2_turns(x, z); + *v = (y + 1.0) * 0.5; +} + +/* Cube: + (x, y, 0) -> (0, x, y) + (x, y, 1) -> (1, x, y) + (0, y, z) -> (2, z, y) + (1, y, z) -> (3, z, y) + (x, 1, z) -> (4, x, z) + (x, 0, z) -> (5, x, z) + Watch out for rounding errors when determining which coordinate is 0 or 1; + see which one is closest to 0 or 1 */ + +static inline void cube_coords(flt x, flt y, flt z, + int * face, flt * u, flt * v) +{ + flt dists[6] = + { fabs(z), fabs(1 - z), fabs(x), fabs(1 - x), fabs(y), fabs(1 - y) }; + flt min; + int f, i; + + f = 0; min = dists[0]; + for (i = 1; i < 6; i++) + if (dists[i] < min) { min = dists[i]; f = i; } + *face = f; + switch (f) { + case 0: case 1: + *u = x; *v = y; break; + case 2: case 3: + *u = z; *v = y; break; + case 4: case 5: + *u = x; *v = z; break; + } +} + +/* Cylinder: + (x, 0, z) -> (2, u, v) where x = 2u-1 and z = 2v-1 hence + (x, 1, z) -> (1, x, z) u = (x+1)/2 and v = (z+1)/2 + (x, y, z) -> (0, u, v) where x = sin(360u) v = y z = cos(360u) + hence u = atan2_turns(x, z) and v = y +*/ + +static inline void cylinder_coords(flt x, flt y, flt z, + int * face, flt * u, flt * v) +{ + flt min, d; + int f; + + min = y * y; f = 0; + d = (1 - y) * (1 - y); + if (d < min) { min = d; f = 1; } + d = fabs(x * x + z * z - 1); + if (d < min) { min = d; f = 2; } + *face = 2 - f; + if (f < 2) { + *u = (x + 1) * 0.5; + *v = (z + 1) * 0.5; + } else { + *u = atan2_turns(x, z); + *v = y; + } +} + +/* Cone: + (x, y, z) -> (0, u, v) where x = v sin 360u + y = v + z = v cos 360u + hence u = atan2_turns(x, z) and v = y + (x, 1, z) -> (1, u, v) where x = 2u-1 and z = 2v-1 + hence u = (x+1)/2 and v = (z+1)/2 +*/ + +static inline void cone_coords(flt x, flt y, flt z, + int * face, flt * u, flt * v) +{ + if ((1 - y) * (1 - y) < fabs(x * x + z * z - y * y)) { + *face = 1; + *u = (x + 1) * 0.5; + *v = (z + 1) * 0.5; + } else { + *face = 0; + *u = atan2_turns(x, z); + *v = y; + } +} + +/* Plane */ + +static inline void plane_coords(flt x, flt y, flt z, + int * face, flt * u, flt * v) +{ + *face = 0; + *u = x; + *v = z; +} + +/* All together */ + +void surface_coords(struct object * obj, struct point * p, + /*out*/ int * face, + /*out*/ flt * u, + /*out*/ flt * v) +{ + switch (obj->kind) { + case Cone: + cone_coords(p->x, p->y, p->z, face, u, v); break; + case Cube: + cube_coords(p->x, p->y, p->z, face, u, v); break; + case Cylinder: + cylinder_coords(p->x, p->y, p->z, face, u, v); break; + case Plane: + plane_coords(p->x, p->y, p->z, face, u, v); break; + case Sphere: + sphere_coords(p->x, p->y, p->z, face, u, v); break; + default: + assert(0); + } +} + diff --git a/test/raytracer/surface.h b/test/raytracer/surface.h new file mode 100644 index 0000000..ac128b2 --- /dev/null +++ b/test/raytracer/surface.h @@ -0,0 +1,7 @@ +/* Return the texture coordinates for the given point on the given + object. Point in object coords. */ + +void surface_coords(struct object * obj, struct point * p, + /*out*/ int * face, + /*out*/ flt * u, + /*out*/ flt * v); diff --git a/test/raytracer/vector.c b/test/raytracer/vector.c new file mode 100644 index 0000000..226b61a --- /dev/null +++ b/test/raytracer/vector.c @@ -0,0 +1,74 @@ +#include "config.h" +#include "point.h" +#include "vector.h" + +flt dotproduct(struct vector * a, struct vector * b) +{ + return a->dx * b->dx + a->dy * b->dy + a->dz * b->dz; +} + +void between(struct point * p, struct point * q, + /*out*/ struct vector * v) +{ + v->dx = q->x - p->x; + v->dy = q->y - p->y; + v->dz = q->z - p->z; +} + +void opposite(struct vector * v, + /*out*/ struct vector * w) +{ + w->dx = - v->dx; + w->dy = - v->dy; + w->dz = - v->dz; +} + +void point_along(struct point * p, struct vector * v, + flt ac, + /*out*/ struct point * q) +{ + q->x = p->x + v->dx * ac; + q->y = p->y + v->dy * ac; + q->z = p->z + v->dz * ac; +} + +void product(struct vector * a, struct vector * b, + /*out*/ struct vector * v) +{ + v->dx = a->dy * b->dz - a->dz * b->dy; + v->dy = a->dz * b->dx - a->dx * b->dz; + v->dz = a->dx * b->dy - a->dy * b->dx; +} + +flt vlength2(struct vector * a) +{ + return a->dx * a->dx + a->dy * a->dy + a->dz * a->dz; +} + +flt vlength(struct vector * a) +{ + return sqrt(vlength2(a)); +} + +void vscale(struct vector * a, flt s, + /*out*/ struct vector * v) +{ + v->dx = a->dx * s; + v->dy = a->dy * s; + v->dz = a->dz * s; +} + +void vnormalize(struct vector * a, /*out*/ struct vector * v) +{ + vscale(a, 1 / vlength(a), v); +} + +void vsub(struct vector * a, struct vector * b, + /*out*/ struct vector * v) +{ + v->dx = a->dx - b->dx; + v->dy = a->dy - b->dy; + v->dz = a->dz - b->dz; +} + + diff --git a/test/raytracer/vector.h b/test/raytracer/vector.h new file mode 100644 index 0000000..c5e400a --- /dev/null +++ b/test/raytracer/vector.h @@ -0,0 +1,23 @@ +struct vector { + flt dx, dy, dz; +}; + +flt dotproduct(struct vector * a, struct vector * b); +void between(struct point * p, struct point * q, + /*out*/ struct vector * v); +void opposite(struct vector * v, + /*out*/ struct vector * w); +void point_along(struct point * p, struct vector * v, + flt ac, + /*out*/ struct point * q); +void product(struct vector * a, struct vector * b, + /*out*/ struct vector * v); +flt vlength2(struct vector * a); +flt vlength(struct vector * a); + +void vscale(struct vector * a, flt s, + /*out*/ struct vector * v); +void vnormalize(struct vector * a, /*out*/ struct vector * v); +void vsub(struct vector * a, struct vector * b, + /*out*/ struct vector * v); + -- cgit v1.2.3