[comp.graphics] ray tracer source

ph@pixar.UUCP (Paul Heckbert) (07/12/87)

The following is C source to a ray tracer written by Masataka Ohta of Japan
which he mailed me in connection with my minimal ray tracer contest.
It was too good (and hence too long) to be a contender in that contest, but
he asked me to post it to the net.
Masataka also wrote the ray tracer used in SEIBU/SEDIC's "Mandala" and
"Maitreya" animations, which were shown at SIGGRAPH a few years back.

		Paul Heckbert

ps: I made minor enhancements to his code to set resolution at run-time.

==========================================

# to unpack, cut here and run the following shell archive through sh
# contents: README Makefile ray.h alloc.c bg.c dindex.c initdir.c
# input.c intersect.c main.c scan.c shade.c shape.c trace.c test.v
# random.c
#
sed 's/^X//' <<'EOF26628' >README
XThis is a public domain ray tracing program with the following features:
X    * shadows
X    * specular reflection
X    * transparency with refraction
X    * antialiasing
X    * high speed
X
XIt computes ray-object intersections very fast by exploiting ray coherence.
XFor many scenes, the program can compute ray-object intersections
Xin constant time, regardless of the number of objects.
X
XThe technique is described in the paper "Ray Coherence Theorem
Xand Constant Time Ray Tracing Algorithm",
X"Computer Graphics 1987, Proceedings of CG International '87"
XSpringer-Verlag, pp. 303-314.
X
XAny suggestions or beautiful pictures generated
Xusing this program are welcome.
X
X			Masataka Ohta
X
X			Computer Center, Tokyo Institute of Technology,
X			2-12-1, O-okayama, Meguro-ku, Tokyo 152, JAPAN
X
X			mohta%titcce.cc.titech.junet@relay.cs.net
X
X==========================================================
X
XTo test the program:
X
X	"make test.pic" will produce a beautiful(?) image of nine spheres.
X
X	"make test1000.pic" will produce an image of 1,000 spheres.
X
X
XUsage: ray [options] <scene file> <picture file>
X
XOptions:
X
X	-a:
X		do anti-aliasing
X
X	-r res
X		set resolution to res (default 128)
X
X	-t:
X		output timing information
X
X	-s:
X		use classical slow algorithm
X	
X	-o:
X		disable object ordering optimization
X
X	-d:
X		debug
X
X	-b:
X		no shading
X
X	
XFile format: 
X	
X	scene file: (an ascii file, see test.v or "make test1000.v" for example)
X		f <field of view>
X		l <x coordinate of the first light source>
X		  <y coordinate of the first light source>
X		  <z coordinate of the first light source>
X		l <x coordinate of the second light source>
X		  <y coordinate of the second light source>
X		  <z coordinate of the second light source>
X		.
X		.
X		.
X		o <shape number (integer) of the first object>
X		  <shade number (integer) of the first object>
X		  <first shape parameter of the first object>
X		   .
X		   .
X		   .
X		  <first shade parameter of the first object>
X		   .
X		   .
X		   .
X		o <shape number (integer) of the second object>
X		  <shade number (integer) of the second object>
X		.
X		.
X		.
X		e
X
X	picture file: (a binary file, make test.pic for an example)
X		<long int xres>
X		<long int yres>
X		<pixel array of yres*xres*3 bytes>
X	pixel array has the form:
X		struct {unsigned char r, g, b;} array[yres][xres];
X	screen space y=0 is at the top
X
XCurrently supported light source types: (should be programmable)
X
X	point light source with constant intensity of (1,1,1)
X	regardless to the distance to the light source
X
XCurrently supported shapes: (may be changed by modifying shape.c)
X
X	number:			0
X	shape:			sphere
X	number of parameters:	4
X		first parameter:	x coordinate value of the center
X		second parameter:	y coordinate value of the center
X		third parameter:	z coordinate value of the center
X		forth parameter:	radius
X
XCurrently supported shades: (may be changed by modifying shade.c)
X
X	number:			0
X	shade:			lambert with shadows
X	number of parameters:	6
X		first to third parameters:	RGB value of ambient
X		4th to 6th parameters:		RGB value of diffuse
X
X	number:			1
X	shade:			phong type specular with shadows
X	number of parameters:	10
X		first to third parameters:	RGB value of ambient
X		4th to 6th parameters:		RGB value of diffuse
X		7th to 9th parameters:		RGB value of specular
X		10th parameter:			specular width factor
X
X	number:			2
X	shade:			mirror (i.e. reflective ) with phong
X				type specular with shadows
X	number of parameters:	13
X		first to third parameters:	RGB value of ambient
X		4th to 6th parameters:		RGB value of diffuse
X		7th to 9th parameters:		RGB value of specular
X		10th parameter:			specular width factor
X		11th to 13th parameters:	RGB ratio of reflection
X
X	number:			3
X	shade:			glass (i.e. reflective and refractive)
X				with phong type specular with shadows
X	number of parameters:	17
X		first to third parameters:	RGB value of ambient
X		4th to 6th parameters:		RGB value of diffuse
X		7th to 9th parameters:		RGB value of specular
X		10th parameter:			specular width factor
X		11th to 13th parameters:	RGB ratio of reflection
X		14th to 16th parameters:	RGB ratio of refraction
X		17th parameter:			refractive index
X
XViewing parameters:
X
X	eye:			fixed at (0,0,0)
X	eye direction:		fixed to (0,0,1)
X	field of view:		variable from 0 to 180 degree (with 'f' option)
X
XResolution:
X
X	set with 'r' option, defaults to 128
X	picture is always square
X
XLimitations:
X
X	coordinate transformation should be supported
X
X	light source type should be programmable
X
X	more complex shapes should be supported
X
X	more sophisticated syntax should be used (eg. comments, symbolic names)
X
X	dynamic loading of light source, shape and shade functions
X	should be possible
EOF26628
sed 's/^X//' <<'EOF26629' >Makefile
XSRC=alloc.c bg.c dindex.c initdir.c input.c intersect.c main.c\
X	scan.c shade.c shape.c trace.c
X
XOBJ=alloc.o bg.o dindex.o initdir.o input.o intersect.o main.o\
X	scan.o shade.o shape.o trace.o
X
XHDR=ray.h
X
XCFLAGS=-O
X
Xray:	$(OBJ)
X	cc $(CFLAGS) $(OBJ) -o ray -lm
X
Xray.lint:	$(HDR) $(SRC)
X	lint -b $(SRC) -lm >ray.lint
X
X$(OBJ):	ray.h
X
Xprint:
X	print Makefile $(HDR) $(SRC)
X
Xray.shar: README Makefile $(HDR) $(SRC) test.v random.c
X	shar README Makefile $(HDR) $(SRC) test.v random.c >ray.shar
X
Xclean:
X	rm -f $(OBJ) ray random test.pic test1000.v test1000.pic lint
X
Xrandom: random.c
X	cc -O random.c -o random -lm
X
Xtest.pic: ray
X	ray test.v test.pic
X
Xtest1000.v: random
X	echo f 90 >test1000.v
X	echo l 1 1 -3 >>test1000.v
X	echo l -4 0 -1 >>test1000.v
X	random 1000 0.07 >>test1000.v
X	echo e >>test1000.v
X
Xtest1000.pic: ray test1000.v
X	ray test1000.v test1000.pic
EOF26629
sed 's/^X//' <<'EOF26630' >ray.h
X#define	MAXOBJECTS	1010
X#define	MAXLEVEL	4
X
X#define	MINT	1e-6
X
Xstruct point {double x,y,z;};
X
Xstruct vector {double x,y,z;};
X
Xstruct circle {
X	struct vector o;
X	double c,s;
X};
X
Xstruct ray {
X	int obj;		/* origin index */
X	struct point a;		/* origin coordinate */
X	struct vector l;	/* direction */
X};
X
Xstruct color {double r,g,b};
X
Xstruct intersect {
X	int obj;
X	double t;
X	struct vector n;
X} intersect();
X
Xstruct object {
X	int flags;
X	struct point center;
X	double radius;
X	struct intersect (*shape)();
X	double *spparams;
X	struct color (*shade)();
X	double *sdparams;
X}
X	objects[MAXOBJECTS];
X
X#define	EYE		1		/* object is eye */
X#define	LIGHT		2		/* object is light source */
X#define	REFLECT		4		/* object is reflective */
X#define	REFRACT		8		/* object is refractive */
X#define	RAYORIGIN	(EYE|LIGHT|REFLECT|REFRACT)
X#define	SELF		16		/* rays from an object may */
X					/* intersect the object again */
X
Xint res;
Xint maxobj,lights;
X
Xstruct objlist2 {
X	int obj;
X	struct objlist *next;
X};
X
Xstruct objlist {
X	int obj;
X	struct objlist *next;
X	float t;
X}
X	***candidates;
X
Xstruct dirlist {
X	int dir;
X	struct dirlist *next;
X};
X
Xstruct vector nilvect,normalize(),reflect(),refract();
Xstruct color trace(),bgcolor();
X
Xint fflag,dflag,aflag,tflag,sflag,oflag;
X
Xint alsize;
X
X/* This should be modified to be word alignment of alloc macro
X *	0x1	for 2 byte alignment
X *	0x3	for 4 byte alignment
X */
X#define	ALMASK	0x1
X
Xchar *myalloc(),*freep,*endp;
Xchar *malloc();
X
X#define alloc(type)	((alsize=(sizeof(type)+ALMASK)&~ALMASK),\
X			 (type *) (endp>freep+alsize\
X				 ?(freep+=alsize)-alsize\
X				 :myalloc(alsize)))
X
X#define nalloc(type,n)	((alsize=(sizeof(type)*(n)+ALMASK)&~ALMASK),\
X			 (type *) (endp>freep+alsize\
X				 ?(freep+=alsize)-alsize\
X				 :myalloc(alsize)))
X
Xint raycount,shadowcount;
X
Xint NRECTS,DIRECTIONS;
X
Xstruct shapetab {
X	void (*shapeinitfunc)();
X	struct intersect (*shapefunc)();
X	int nparams;
X}
X	*shapetab;
X
Xint nshapetab;
X
Xstruct shadetab {
X	void (*shadeinitfunc)();
X	struct color (*shadefunc)();
X	int nparams;
X}
X	*shadetab;
X
Xint nshadetab;
X
Xdouble fovf;
EOF26630
sed 's/^X//' <<'EOF26631' >alloc.c
X#include "ray.h"
X
Xchar *freep=(char *)0,*endp=(char *)0,*malloc();
X
X#define ALLOCUNIT 16384
X
Xchar *myalloc(s)
Xint s;
X{int size;
Xchar *a;
X	if(s>ALLOCUNIT)
X		size=s;
X	else
X		size=ALLOCUNIT;
X	a=malloc((unsigned)size);
X	if(a==0)
X	{	perror("malloc");
X		exit(1);
X	}
X	freep=a+s;
X	endp=a+size;
X	return a;
X}
EOF26631
sed 's/^X//' <<'EOF26632' >bg.c
X#include "ray.h"
X
X/*ARGSUSED*/
Xstruct color bgcolor(r)
Xstruct ray r;
X{struct color c;
X	if(r.l.y>0)
X		c.r=c.g=c.b=0.2;
X	else
X		c.r=c.g=c.b=0.2-0.8*r.l.y;
X	return c;
X}
EOF26632
sed 's/^X//' <<'EOF26633' >dindex.c
X#include "ray.h"
X
Xdindex(v)
Xstruct vector v;
X{double x,y,z,ax,ay,az;
X	x=v.x;
X	y=v.y;
X	z=v.z;
X	ax=(x>0)?x:-x;
X	ay=(y>0)?y:-y;
X	az=(z>0)?z:-z;
X	if(ax>=ay&&ax>=az)
X		if(x>0)
X			return index2(y,z,ax);
X		else
X			return index2(y,z,ax)+(NRECTS*NRECTS);
X	else if(ay>=ax&&ay>=az)
X		if(y>0)
X			return index2(z,x,ay)+(2*NRECTS*NRECTS);
X		else
X			return index2(z,x,ay)+(3*NRECTS*NRECTS);
X	else
X		if(z>0)
X			return index2(x,y,az)+(4*NRECTS*NRECTS);
X		else
X			return index2(x,y,az)+(5*NRECTS*NRECTS);
X}
X
Xindex2(x,y,r)
Xdouble r,x,y;
X{register int m,n;
X	m=(r+x)/(r*2)*NRECTS;
X	if(m==NRECTS)
X		m=NRECTS-1;
X	n=(r+y)/(r*2)*NRECTS;
X	if(n==NRECTS)
X		n=NRECTS-1;
X	return n*NRECTS+m;
X}
EOF26633
sed 's/^X//' <<'EOF26634' >initdir.c
X#include <math.h>
X#include <stdio.h>
X
X#include "ray.h"
X
Xstruct circle *dvl;
X
Xstruct dirlist **neighbor;
X
Xint *markdir;
X
Xinitdir()
X{register int o0,o1,d;
Xregister double l;
Xstruct point c0,c1;
Xregister double r0,r1;
Xstruct circle cir;
Xregister struct objlist *ol;
Xint origins;
X	origins=0;
X	for(o0=0;o0<maxobj;o0++)
X	{	if((objects[o0].flags&RAYORIGIN))
X			origins++;
X	}
X	NRECTS=40;
X	if(NRECTS*sqrt((double)origins)>100)
X	{	if(origins)
X			NRECTS=100/sqrt((double)origins);
X		else
X			NRECTS=100;
X	}
X	if(NRECTS<4)
X		NRECTS=4;
X	DIRECTIONS=6*NRECTS*NRECTS;
X	dvl=nalloc(struct circle,DIRECTIONS);
X	neighbor=nalloc(struct dirlist *,DIRECTIONS);
X	markdir=nalloc(int,DIRECTIONS);
X	candidates=nalloc(struct objlist **,maxobj);
X	for(o0=0;o0<maxobj;o0++)
X	{	if(!(objects[o0].flags&RAYORIGIN))
X			continue;;
X		candidates[o0]=nalloc(struct objlist *,DIRECTIONS);
X	}
X	initdvl();
X	initneighbor();
X	for(d=0;d<DIRECTIONS;d++)
X		markdir[d]=0;
X	for(o0=0;o0<maxobj;o0++)
X	{	if(!(objects[o0].flags&RAYORIGIN))
X			continue;;
X		for(d=0;d<DIRECTIONS;d++)
X			candidates[o0][d]=(struct objlist *)0;
X		for(o1=lights+1;o1<maxobj;o1++)
X		{	if(o0==o1)
X				continue;
X			if(dflag)
X				fprintf(stderr,"object %d to %d\n",o0,o1);
X			c0=objects[o0].center;
X			r0=objects[o0].radius;
X			c1=objects[o1].center;
X			r1=objects[o1].radius;
X			cir.o.x=c1.x-c0.x;
X			cir.o.y=c1.y-c0.y;
X			cir.o.z=c1.z-c0.z;
X			l=sqrt(cir.o.x*cir.o.x+cir.o.y*cir.o.y+cir.o.z*cir.o.z);
X			cir.o.x/=l;
X			cir.o.y/=l;
X			cir.o.z/=l;
X			if(r0+r1>=l)
X			{	for(d=0;d<DIRECTIONS;d++)
X				{	if(oflag)
X						ol=alloc(struct objlist);
X					else
X						ol=(struct objlist *)
X							alloc(struct objlist2);
X					ol->obj=o1;
X					ol->next=candidates[o0][d];
X					if(oflag)
X						ol->t=0;
X					candidates[o0][d]=ol;
X				}
X			}
X			else
X			{	cir.s=(r0+r1)/l;
X				cir.c=sqrt(1-cir.s*cir.s);
X				d=dindex(cir.o);
X				if(oflag)
X					ttraverse(d,&cir,o1,
X						candidates[o0],l-r0-r1);
X				else
X					traverse(d,&cir,o1,candidates[o0]);
X				tclear(d);
X			}
X		}
X		if(oflag)
X			for(d=0;d<DIRECTIONS;d++)
X				sortlist(candidates[o0][d]);
X	}
X	if(tflag)
X	{int c;
X		c=0;
X		for(d=0;d<DIRECTIONS;d++)
X			for(ol=candidates[0][d];ol;ol=ol->next)
X				c++;
X		fprintf(stderr,"Image complexity %g\n",((double)c)/DIRECTIONS);
X		for(o0=1;o0<lights+1;o0++)
X		{	c=0;
X			for(d=0;d<DIRECTIONS;d++)
X				for(ol=candidates[o0][d];ol;ol=ol->next)
X					c++;
X			fprintf(stderr,"Shadow complexity for light %d %g\n",
X				o0,((double)c)/DIRECTIONS);
X		}
X	}
X}
X
Xtraverse(d,cir,o,ol)
Xregister int d;
Xregister struct circle *cir;
Xregister int o;
Xregister struct objlist **ol;
X{register struct dirlist *dl;
Xregister struct objlist *ol2;
X	markdir[d]=1;
X	if(!overlap(&dvl[d],cir))
X		return;
X	ol2=(struct objlist *)alloc(struct objlist2);
X	ol2->obj=o;
X	ol2->next=ol[d];
X	ol[d]=ol2;
X	for(dl=neighbor[d];dl;dl=dl->next)
X		if(!markdir[dl->dir])
X			traverse(dl->dir,cir,o,ol);
X}
X
Xttraverse(d,cir,o,ol,t)
Xregister int d;
Xregister struct circle *cir;
Xregister int o;
Xregister struct objlist **ol;
Xregister double t;
X{register struct dirlist *dl;
Xregister struct objlist *ol2;
X	markdir[d]=1;
X	if(!overlap(&dvl[d],cir))
X		return;
X	ol2=alloc(struct objlist);
X	ol2->obj=o;
X	ol2->next=ol[d];
X	ol2->t=t;
X	ol[d]=ol2;
X	for(dl=neighbor[d];dl;dl=dl->next)
X		if(!markdir[dl->dir])
X			ttraverse(dl->dir,cir,o,ol,t);
X}
X
Xtclear(d)
Xregister int d;
X{register struct dirlist *dl;
X	markdir[d]=0;
X	for(dl=neighbor[d];dl;dl=dl->next)
X		if(markdir[dl->dir])
X			tclear(dl->dir);
X}
X
Xoverlap(cir0,cir1)
Xregister struct circle *cir0,*cir1;
X{	return
X		cir0->o.x*cir1->o.x+cir0->o.y*cir1->o.y+cir0->o.z*cir1->o.z
X		>
X		cir0->c*cir1->c-cir0->s*cir1->s;
X}
X
Xstruct circle initdvl2(x,y,z,dx0,dy0,dz0,dx1,dy1,dz1)
Xdouble x,y,z,dx0,dy0,dz0,dx1,dy1,dz1;
X{int i;
Xdouble l,c,minc;
Xstruct vector v[4];
Xstruct circle cir;
X	for(i=0;i<4;i++)
X	{	v[i].x=x;
X		v[i].y=y;
X		v[i].z=z;
X	}
X	v[1].x+=dx0;
X	v[1].y+=dy0;
X	v[1].z+=dz0;
X	v[2].x+=dx0+dx1;
X	v[2].y+=dy0+dy1;
X	v[2].z+=dz0+dz1;
X	v[3].x+=dx1;
X	v[3].y+=dy1;
X	v[3].z+=dz1;
X	cir.o.x=cir.o.y=cir.o.z=0;
X	for(i=0;i<4;i++)
X	{	cir.o.x+=v[i].x/4;
X		cir.o.y+=v[i].y/4;
X		cir.o.z+=v[i].z/4;
X		l=1/sqrt(v[i].x*v[i].x+v[i].y*v[i].y+v[i].z*v[i].z);
X		v[i].x*=l;
X		v[i].y*=l;
X		v[i].z*=l;
X	}
X	l=1/sqrt(cir.o.x*cir.o.x+cir.o.y*cir.o.y+cir.o.z*cir.o.z);
X	cir.o.x*=l;
X	cir.o.y*=l;
X	cir.o.z*=l;
X	minc=1;
X	for(i=0;i<4;i++)
X	{	c=cir.o.x*v[i].x+cir.o.y*v[i].y+cir.o.z*v[i].z;
X		if(c<minc)
X			minc=c;
X	}
X	cir.c=minc*0.999999;
X	cir.s=sqrt(1-cir.c*cir.c);
X	return cir;
X}
X
Xinitdvl()
X{int i,j,d;
Xdouble x,y,z,dl;
X	dl=2.0/NRECTS;
X	d=0;
X	x=1;
X	z= -1;
X	for(i=0;i<NRECTS;i++)
X	{	y= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,0.0,0.0,dl,0.0,dl,0.0);
X			y+=dl;
X		}
X		z+=dl;
X	}
X	x= -1;
X	z= -1;
X	for(i=0;i<NRECTS;i++)
X	{	y= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,0.0,dl,0.0,0.0,0.0,dl);
X			y+=dl;
X		}
X		z+=dl;
X	}
X	y=1;
X	x= -1;
X	for(i=0;i<NRECTS;i++)
X	{	z= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,dl,0.0,0.0,0.0,0.0,dl);
X			z+=dl;
X		}
X		x+=dl;
X	}
X	y= -1;
X	x= -1;
X	for(i=0;i<NRECTS;i++)
X	{	z= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,0.0,0.0,dl,dl,0.0,0.0);
X			z+=dl;
X		}
X		x+=dl;
X	}
X	z=1;
X	y= -1;
X	for(i=0;i<NRECTS;i++)
X	{	x= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,0.0,dl,0.0,dl,0.0,0.0);
X			x+=dl;
X		}
X		y+=dl;
X	}
X	z= -1;
X	y= -1;
X	for(i=0;i<NRECTS;i++)
X	{	x= -1;
X		for(j=0;j<NRECTS;j++)
X		{	dvl[d++]=initdvl2(x,y,z,dl,0.0,0.0,0.0,dl,0.0);
X			x+=dl;
X		}
X		y+=dl;
X	}
X}
X
Xstruct dirlist *consdir(m,dl)
Xint m;
Xstruct dirlist *dl;
X{struct dirlist *dl2;
X	dl2=alloc(struct dirlist);
X	dl2->dir=m;
X	dl2->next=dl;
X	return dl2;
X}
X
Xinitneighbor()
X{int i;
Xdouble d;
X	for(i=0;i<DIRECTIONS;i++)
X		neighbor[i]=(struct dirlist *)0;
X	d=2.0/NRECTS;
X	initneighbor2(1.0,-1.0,-1.0,0.0,d,0.0,0.0,0.0,d);
X	initneighbor2(-1.0,-1.0,-1.0,0.0,d,0.0,0.0,0.0,d);
X	initneighbor2(-1.0,1.0,-1.0,d,0.0,0.0,0.0,0.0,d);
X	initneighbor2(-1.0,-1.0,-1.0,d,0.0,0.0,0.0,0.0,d);
X	initneighbor2(-1.0,-1.0,1.0,d,0.0,0.0,0.0,d,0.0);
X	initneighbor2(-1.0,-1.0,-1.0,d,0.0,0.0,0.0,d,0.0);
X}
X
Xinitneighbor2(x,y,z,dx0,dy0,dz0,dx1,dy1,dz1)
Xdouble x,y,z,dx0,dy0,dz0,dx1,dy1,dz1;
X{int i,j,d;
Xdouble x0,y0,z0,x1,y1,z1;
X	x0=x+dx0/2+dx1/2;
X	y0=y+dy0/2+dy1/2;
X	z0=z+dz0/2+dz1/2;
X	for(i=0;i<NRECTS;i++)
X	{	for(j=0;j<NRECTS;j++)
X		{	x1=x0+i*dx0+j*dx1;
X			y1=y0+i*dy0+j*dy1;
X			z1=z0+i*dz0+j*dz1;
X			d=dindex2(x1,y1,z1);
X			neighbor[d]=consdir(dindex2(x1+dx0,y1+dy0,z1+dz0),
X				neighbor[d]);
X			neighbor[d]=consdir(dindex2(x1-dx0,y1-dy0,z1-dz0),
X				neighbor[d]);
X			neighbor[d]=consdir(dindex2(x1+dx1,y1+dy1,z1+dz1),
X				neighbor[d]);
X			neighbor[d]=consdir(dindex2(x1-dx1,y1-dy1,z1-dz1),
X				neighbor[d]);
X		}
X	}
X}
X
Xdindex2(x,y,z)
Xdouble x,y,z;
X{struct vector v;
X	v.x=x;
X	v.y=y;
X	v.z=z;
X	return dindex(v);
X}
X
Xprint2(v)
Xstruct vector *v;
X{	if(v->z==0)
X	{	printf("\tz=0\n");
X		return;
X	}
X	printf("\t%g %g\n",v->x/v->z,v->y/v->z);
X}
Xsortlist(ol)
Xregister struct objlist *ol;
X{register struct objlist *ol2,*ol3;
Xdouble t;
Xint obj;
X	if(!ol)
X		return;
X	for(;ol->next;ol=ol->next)
X	{	ol3=ol;
X		for(ol2=ol->next;ol2;ol2=ol2->next)
X			if(ol2->t<ol3->t)
X				ol3=ol2;
X		if(ol!=ol3)
X		{	t=ol->t;
X			ol->t=ol3->t;
X			ol3->t=t;
X			obj=ol->obj;
X			ol->obj=ol3->obj;
X			ol3->obj=obj;
X		}
X	}
X}
EOF26634
sed 's/^X//' <<'EOF26635' >input.c
X#include <math.h>
X#include <stdio.h>
X
X#include "ray.h"
X
Xinput(s)
Xchar *s;
X{struct object *o;
XFILE *f;
Xint i,sp,sd;
Xchar c;
Xdouble fov;
X	if((f=fopen(s,"r"))==0)
X	{	perror(s);
X		exit(1);
X	}
X	o=objects;
X	if(fscanf(f," f %lf ",&fov)!=1||fov<=0||fov>=180)
X	{	fprintf(stderr,"field of view format error\n");
X		exit(1);
X	}
X	fovf=tan(fov*3.14159/360)/sqrt(2.0);
X	o->flags=EYE;
X	o->center.x=o->center.y=o->center.z=0;
X	o->radius=0;
X	o++;
X	lights=0;
X	while(fscanf(f," l %lf %lf %lf ",&o->center.x,&o->center.y,
X		&o->center.z)==3)
X	{	o->flags=LIGHT;
X		o->radius=0;
X		o++;
X		lights++;
X	}
X	maxobj=lights+1;
X	while(fscanf(f," o %d %d ",&sp,&sd)==2)
X	{	if(maxobj==MAXOBJECTS)
X		{	fprintf("too much objects\n");
X			exit(1);
X		}
X		if(sp<0||sp>=nshapetab)
X		{	fprintf(stderr,"shape number range error\n");
X			exit(1);
X		}
X		if(sd<0||sd>=nshadetab)
X		{	fprintf(stderr,"shade number range error\n");
X			exit(1);
X		}
X		o= &objects[maxobj];
X		o->flags=0;
X
X		o->shape=shapetab[sp].shapefunc;
X		o->spparams=nalloc(double,shapetab[sp].nparams);
X		for(i=0;i<shapetab[sp].nparams;i++)
X		{	if(fscanf(f," %lf ",&o->spparams[i])!=1)
X			{	fprintf(stderr,"shape parameter error\n");
X				exit(1);
X			}
X		}
X		(*shapetab[sp].shapeinitfunc)(o);
X
X		o->shade=shadetab[sd].shadefunc;
X		o->sdparams=nalloc(double,shadetab[sd].nparams);
X		for(i=0;i<shadetab[sd].nparams;i++)
X		{	if(fscanf(f," %lf ",&o->sdparams[i])!=1)
X			{	fprintf(stderr,"shade parameter error\n");
X				exit(1);
X			}
X		}
X		(*shadetab[sd].shadeinitfunc)(o);
X
X		maxobj++;
X	}
X	if(fscanf(f," %c",&c)!=1||c!='e')
X	{	fprintf(stderr,"end command not found\n");
X		exit(1);
X	}
X	(void) fclose(f);
X}
EOF26635
sed 's/^X//' <<'EOF26636' >intersect.c
X#include "ray.h"
X
Xstruct intersect intersect(r)
Xstruct ray r;
X{struct intersect i,imin;
Xregister struct objlist *ol;
Xregister int o;
X	imin.obj=0;
X	if(fflag)
X	{	if(objects[r.obj].flags&SELF)
X		{	i=(*objects[r.obj].shape)(r,r.obj);
X			if(i.obj)
X				imin=i;
X		}
X		for(ol=candidates[r.obj][dindex(r.l)];ol;ol=ol->next)
X		{	if(oflag&&imin.obj&&imin.t<ol->t)
X				break;
X			i=(*objects[ol->obj].shape)(r,ol->obj);
X			if(i.obj&&(!imin.obj||i.t<imin.t))
X				imin=i;
X		}
X	}
X	else
X		for(o=lights+1;o<maxobj;o++)
X		{	i=(*objects[o].shape)(r,o);
X			if(i.obj&&(!imin.obj||i.t<imin.t))
X				imin=i;
X		}
X	return imin;
X}
EOF26636
sed 's/^X//' <<'EOF26637' >main.c
X#include <sys/time.h>
X#include <sys/resource.h>
X#include <stdio.h>
X
X#include "ray.h"
X
Xmain(ac,av)
Xint ac;
Xchar **av;
X{int g;
Xstruct rusage ru;
Xstruct timeval tu,tu2;
X	raycount=0;
X	shadowcount=0;
X	fflag=1;
X	dflag=0;
X	aflag=0;
X	tflag=0;
X	sflag=1;
X	oflag=1;
X	while(av[1]&&*av[1]=='-')
X	{	switch(av[1][1])
X		{case 's':		/* classical slow algorithm */
X			fflag=0;
X			break;
X		case 'd':		/* debug flag */
X			dflag=1;
X			break;
X		case 'a':		/* do anti-aliasing */
X			aflag=1;
X			break;
X		case 't':		/* produce timing information */
X			tflag=1;
X			break;
X		case 'b':		/* no shading */
X			sflag=0;
X			break;
X		case 'o':		/* don't order candidates */
X			oflag=0;
X			break;
X		case 'r':		/* set image resolution */
X			res = atoi(av[2]);
X			ac--;
X			av++;
X			break;
X		default:
X			usage();
X		}
X		ac--;
X		av++;
X	}
X	if(ac!=3)
X		usage();
X	initshape();
X	initshade();
X	input(av[1]);
X	if(fflag)
X	{	(void) getrusage(0,&ru);
X		tu=ru.ru_utime;
X		initdir();
X		if(tflag)
X			fprintf(stderr,"Directions=%d Objects=%d\n",
X				DIRECTIONS,maxobj);
X		(void) getrusage(0,&ru);
X		tu2=ru.ru_utime;
X		if(tflag)
X		{	fprintf(stderr,"Initializeation time ");
X			printtime(tu,tu2);
X		}
X	}
X	else
X	{	if(tflag)
X			fprintf(stderr,"Objects=%d\n",maxobj);
X	}
X	if((g=creat(av[2],0666))<0)
X	{	perror(av[2]);
X		exit(3);
X	}
X	(void) getrusage(0,&ru);
X	tu=ru.ru_utime;
X	if(aflag)
X		ascan(g);
X	else
X		scan(g);
X	(void) getrusage(0,&ru);
X	tu2=ru.ru_utime;
X	if(tflag)
X	{	fprintf(stderr,"   Computation time ");
X		printtime(tu,tu2);
X		fprintf(stderr,"\n%d rays traced",raycount);
X		fprintf(stderr,"   %d shadows tested\n",shadowcount);
X	}
X	(void) close(g);
X	exit(0);
X}
X
Xusage()
X{	fprintf(stderr,
X		"Usage: ray [-sdatbo] [-r res] scenefile picfile\n");
X	exit(1);
X}
X
Xprinttime(t0,t1)
Xstruct timeval t0,t1;
X{int sec,usec;
X	usec=t1.tv_usec-t0.tv_usec;
X	sec=t1.tv_sec-t0.tv_sec;
X	if(usec<0)
X	{	usec+=1000000;
X		sec--;
X	}
X	fprintf(stderr,"%4d.%03d",sec,usec/1000);
X}
EOF26637
sed 's/^X//' <<'EOF26638' >scan.c
X#include <stdio.h>
X#include <math.h>
X
X#include "ray.h"
X
Xint res = 128, xres, yres;
X#define MAXRES 2048
X
Xtypedef struct {int xres, yres;} pichead;
Xtypedef unsigned char icolor[3];
X
Xchar aaflag[5][4*MAXRES+1];
Xicolor aabuf[5][4*MAXRES+1];
X
Xpixel(ic,x,y)
Xicolor ic;
Xregister double x,y;
X{register double t;
Xstruct ray r;
Xstruct color c;
Xregister double m;
X	x*=fovf;
X	y*=fovf;
X	r.obj=0;
X	r.a.x=r.a.y=r.a.z=0;
X	t=1/sqrt(x*x+y*y+1);
X	r.l.x=x*t;
X	r.l.y=y*t;
X	r.l.z=t;
X	c=trace(0,r);
X	m=c.r;
X	if(m<c.g)
X		m=c.g;
X	if(m<c.b)
X		m=c.b;
X	if(m>1.0)
X	{	c.r/=m;
X		c.g/=m;
X		c.b/=m;
X	}
X	ic[0]=c.r*255;
X	ic[1]=c.g*255;
X	ic[2]=c.b*255;
X}
X
Xaapixel2(ic,i,j)
Xicolor ic;
Xregister int i,j;
X{register double x,y;
X	x=((double)(j-2*xres))/(2*xres);
X	y=((double)(2*yres-i))/(2*yres);
X	pixel(ic,x,y);
X}
X
Xaapixel(ic,i,j,ai,s)
Xicolor ic;
Xregister int i,j,ai,s;
X{register int k;
Xicolor c,c00,c01,c10,c11;
X	if(!aaflag[ai][j])
X	{	aapixel2(aabuf[ai][j],i,j);
X		aaflag[ai][j]=1;
X	}
X	if(!aaflag[ai+s][j])
X	{	aapixel2(aabuf[ai+s][j],i+s,j);
X		aaflag[ai+s][j]=1;
X	}
X	if(!aaflag[ai][j+s])
X	{	aapixel2(aabuf[ai][j+s],i,j+s);
X		aaflag[ai][j+s]=1;
X	}
X	if(!aaflag[ai+s][j+s])
X	{	aapixel2(aabuf[ai+s][j+s],i+s,j+s);
X		aaflag[ai+s][j+s]=1;
X	}
X	for(k=0;k<3;k++)
X	{	c00[k]=aabuf[ai][j][k];
X		c10[k]=aabuf[ai+s][j][k];
X		c01[k]=aabuf[ai][j+s][k];
X		c11[k]=aabuf[ai+s][j+s][k];
X	}
X	c[0]=(c00[0]+c01[0]+c10[0]+c11[0])/4;
X	c[1]=(c00[1]+c01[1]+c10[1]+c11[1])/4;
X	c[2]=(c00[2]+c01[2]+c10[2]+c11[2])/4;
X	if(s==1||nearc(c,c00)&&nearc(c,c01)&&nearc(c,c10)&&nearc(c,c11))
X	{	ic[0]=c[0];
X		ic[1]=c[1];
X		ic[2]=c[2];
X		return;
X	}
X	s/=2;
X	aapixel(c00,i,j,ai,s);
X	aapixel(c01,i+s,j,ai+s,s);
X	aapixel(c10,i,j+s,ai,s);
X	aapixel(c11,i+s,j+s,ai+s,s);
X	ic[0]=(c00[0]+c01[0]+c10[0]+c11[0])/4;
X	ic[1]=(c00[1]+c01[1]+c10[1]+c11[1])/4;
X	ic[2]=(c00[2]+c01[2]+c10[2]+c11[2])/4;
X}
X
Xscan(f)
Xregister int f;
X{register int i,j;
Xregister double x,y;
Xicolor buf[MAXRES];
X	scaninit(f);
X	for(i=0;i<yres;i++)
X	{	for(j=0;j<xres;j++)
X		{	x=(2*j-xres+0.5)/xres;
X			y=(yres-2*i-0.5)/yres;
X			pixel(buf[j],x,y);
X		}
X		(void) write(f,(char *)buf,xres*sizeof(icolor));
X	}
X}
X
Xascan(f)
Xregister int f;
X{icolor buf[MAXRES];
Xregister double x,y;
Xregister int i,j,ai;
X	scaninit(f);
X	for(ai=0;ai<5;ai++)
X		for(j=0;j<=4*xres;j++)
X			aaflag[ai][j]=0;
X	for(j=0;j<=4*xres;j+=4)
X	{	x=((double)(j-2*xres))/(2*xres);
X		y=1;
X		pixel(aabuf[0][j],x,y);
X		aaflag[0][j]=1;
X	}
X	for(i=0;i<4*yres;i+=4)
X	{	for(j=0;j<=4*xres;j+=4)
X		{	x=((double)(j-2*xres))/(2*xres);
X			y=((double)(2*yres-(i+4)))/(2*yres);
X			pixel(aabuf[4][j],x,y);
X			aaflag[4][j]=1;
X		}
X		for(j=0;j<4*xres;j+=4)
X			aapixel(buf[j/4],i,j,0,4);
X		(void) write(f,(char *)buf,xres*sizeof(icolor));
X		for(j=0;j<=4*xres;j++)
X		{	aabuf[0][j][0]=aabuf[4][j][0];
X			aabuf[0][j][1]=aabuf[4][j][1];
X			aabuf[0][j][2]=aabuf[4][j][2];
X			aaflag[0][j]=aaflag[4][j];
X		}
X		for(ai=1;ai<5;ai++)
X		{	for(j=0;j<=4*xres;j++)
X				aaflag[ai][j]=0;
X		}
X	}
X}
X
Xscaninit(f)
Xint f;
X{pichead head;
X	if (res>MAXRES)
X	{	fprintf(stderr, "res=%d too large, must be less than %d\n",
X			res, MAXRES);
X		exit(1);
X	}
X	head.xres = xres = res;
X	head.yres = yres = res;
X	(void) write(f,(char *)&head,sizeof head);
X}
X
Xnearc(c0,c1)
Xicolor c0,c1;
X{register int d;
X	d=abs(c0[0]-c1[0])+abs(c0[1]-c1[1])+abs(c0[2]-c1[2]);
X	return d<15;
X}
EOF26638
sed 's/^X//' <<'EOF26639' >shade.c
X#include <math.h>
X
X#include "ray.h"
X
Xvoid initdiffuse(),initspecular(),initmirror(),initglass();
Xstruct color diffuse(),specular(),mirror(),glass();
X
Xdouble phong();
X
Xinitshade()
X{	nshadetab=4;
X	shadetab=nalloc(struct shadetab,nshadetab);
X	shadetab[0].shadeinitfunc=initdiffuse;
X	shadetab[0].shadefunc=diffuse;
X	shadetab[0].nparams=6;
X	shadetab[1].shadeinitfunc=initspecular;
X	shadetab[1].shadefunc=specular;
X	shadetab[1].nparams=10;
X	shadetab[2].shadeinitfunc=initmirror;
X	shadetab[2].shadefunc=mirror;
X	shadetab[2].nparams=13;
X	shadetab[3].shadeinitfunc=initglass;
X	shadetab[3].shadefunc=glass;
X	shadetab[3].nparams=17;
X}
X
Xshadow(o,r)
Xint o;
Xstruct ray r;
X{struct intersect i;
X	shadowcount++;
X	i=intersect(r);
X	return i.obj!=o;
X}
X
Xdouble phong(n,e,li,g)
Xstruct vector n,e,li;
Xdouble g;
X{struct vector h;
Xdouble l,c,t;
X	h.x=e.x+li.x;
X	h.y=e.y+li.y;
X	h.z=e.z+li.z;
X	l= -(h.x*n.x+h.y*n.y+h.z*n.z);
X	if(l<0)
X		return 0.0;
X	c=l*l/(h.x*h.x+h.y*h.y+h.z*h.z);
X	t=(1-c)/c;
X	if(t>10*g)
X		return 0;
X	return exp(-t/g);
X}
X
Xstruct vector nilvect={0.0,0.0,0.0};
X
Xstruct vector reflect(n,u)
Xstruct vector n,u;
X{struct vector v;
Xdouble p;
X	p=2*(n.x*u.x+n.y*u.y+n.z*u.z);
X	v.x=u.x-p*n.x;
X	v.y=u.y-p*n.y;
X	v.z=u.z-p*n.z;
X	return v;
X}
X
Xstruct vector refract(n,u,ni)
Xstruct vector n,u;
Xdouble ni;
X{struct vector v;
Xdouble p,t;
X	p=n.x*u.x+n.y*u.y+n.z*u.z;
X	if(p<0)
X	{	t=1-(1-p*p)/(ni*ni);
X		if(t<=0)
X			return nilvect;
X		t= -p/ni-sqrt(t);
X	}
X	else
X	{	ni=1/ni;
X		t=1-(1-p*p)/(ni*ni);
X		if(t<=0)
X			return nilvect;
X		t= -p/ni+sqrt(t);
X	}
X	v.x=u.x/ni+t*n.x;
X	v.y=u.y/ni+t*n.y;
X	v.z=u.z/ni+t*n.z;
X	return v;
X}
X
Xvoid initdiffuse(o)
Xstruct object *o;
X{
X#ifdef lint
X	o=o;
X#endif
X}
X
Xvoid initspecular(o)
Xstruct object *o;
X{
X#ifdef lint
X	o=o;
X#endif
X}
X
Xvoid initmirror(o)
Xstruct object *o;
X{	o->flags|=REFLECT;
X}
X
Xvoid initglass(o)
Xstruct object *o;
X{	o->flags|=REFLECT|REFRACT|SELF;
X}
X
Xstruct color diffuse(n,r,i)
Xint n;
Xstruct ray r;
Xstruct intersect i;
X{int j;
Xstruct color c;
Xdouble l,vx,vy,vz;
Xstruct ray r2;
Xdouble *sdparams;
X#ifdef lint
X	n=n;
X#endif
X	sdparams=objects[i.obj].sdparams;
X	c.r=sdparams[0];
X	c.g=sdparams[1];
X	c.b=sdparams[2];
X	for(j=1;j<lights+1;j++)
X	{	r2.a=objects[j].center;
X		vx=r.a.x+i.t*r.l.x-r2.a.x;
X		vy=r.a.y+i.t*r.l.y-r2.a.y;
X		vz=r.a.z+i.t*r.l.z-r2.a.z;
X		l=1/sqrt(vx*vx+vy*vy+vz*vz);
X		vx*=l;
X		vy*=l;
X		vz*=l;
X		r2.obj=j;
X		r2.l.x=vx;
X		r2.l.y=vy;
X		r2.l.z=vz;
X		l= -(i.n.x*vx+i.n.y*vy+i.n.z*vz);
X		if(l>0&&!shadow(i.obj,r2))
X		{	c.r+=l*sdparams[3];
X			c.g+=l*sdparams[4];
X			c.b+=l*sdparams[5];
X		}
X	}
X	return c;
X}
X
Xstruct color specular(n,r,i)
Xint n;
Xstruct ray r;
Xstruct intersect i;
X{int j;
Xstruct color c;
Xdouble l,x,y,z,sp;
Xstruct ray r2;
Xstruct vector e,li;
Xdouble *sdparams;
X#ifdef lint
X	n=n;
X#endif
X	sdparams=objects[i.obj].sdparams;
X	c.r=sdparams[0];
X	c.g=sdparams[1];
X	c.b=sdparams[2];
X	x=r.a.x+i.t*r.l.x;
X	y=r.a.y+i.t*r.l.y;
X	z=r.a.z+i.t*r.l.z;
X	e.x=x-objects[0].center.x;
X	e.y=y-objects[0].center.y;
X	e.z=z-objects[0].center.z;
X	l=1/sqrt(e.x*e.x+e.y*e.y+e.z*e.z);
X	e.x*=l;
X	e.y*=l;
X	e.z*=l;
X	for(j=1;j<lights+1;j++)
X	{	r2.a=objects[j].center;
X		li.x=x-r2.a.x;
X		li.y=y-r2.a.y;
X		li.z=z-r2.a.z;
X		l=1/sqrt(li.x*li.x+li.y*li.y+li.z*li.z);
X		li.x*=l;
X		li.y*=l;
X		li.z*=l;
X		r2.obj=j;
X		r2.l=li;
X		l= -(i.n.x*li.x+i.n.y*li.y+i.n.z*li.z);
X		if(l>0&&!shadow(i.obj,r2))
X		{	c.r+=l*sdparams[3];
X			c.g+=l*sdparams[4];
X			c.b+=l*sdparams[5];
X			sp=phong(i.n,e,li,sdparams[9]);
X			c.r+=sp*sdparams[6];
X			c.g+=sp*sdparams[7];
X			c.b+=sp*sdparams[8];
X		}
X	}
X	return c;
X}
X
Xstruct color mirror(n,r,i)
Xint n;
Xstruct ray r;
Xstruct intersect i;
X{int j;
Xstruct color c,rc;
Xdouble l,x,y,z,sp;
Xstruct ray r2;
Xstruct vector e,li;
Xdouble *sdparams;
X	sdparams=objects[i.obj].sdparams;
X	c.r=sdparams[0];
X	c.g=sdparams[1];
X	c.b=sdparams[2];
X	x=r.a.x+i.t*r.l.x;
X	y=r.a.y+i.t*r.l.y;
X	z=r.a.z+i.t*r.l.z;
X	e.x=x-objects[0].center.x;
X	e.y=y-objects[0].center.y;
X	e.z=z-objects[0].center.z;
X	l=1/sqrt(e.x*e.x+e.y*e.y+e.z*e.z);
X	e.x*=l;
X	e.y*=l;
X	e.z*=l;
X	for(j=1;j<lights+1;j++)
X	{	r2.a=objects[j].center;
X		li.x=x-r2.a.x;
X		li.y=y-r2.a.y;
X		li.z=z-r2.a.z;
X		l=1/sqrt(li.x*li.x+li.y*li.y+li.z*li.z);
X		li.x*=l;
X		li.y*=l;
X		li.z*=l;
X		r2.obj=j;
X		r2.l=li;
X		l= -(i.n.x*li.x+i.n.y*li.y+i.n.z*li.z);
X		if(l>0&&!shadow(i.obj,r2))
X		{	c.r+=l*sdparams[3];
X			c.g+=l*sdparams[4];
X			c.b+=l*sdparams[5];
X			sp=phong(i.n,e,li,sdparams[9]);
X			c.r+=sp*sdparams[6];
X			c.g+=sp*sdparams[7];
X			c.b+=sp*sdparams[8];
X		}
X	}
X	r2.obj=i.obj;
X	r2.a.x=x;
X	r2.a.y=y;
X	r2.a.z=z;
X	r2.l=reflect(i.n,r.l);
X	rc=trace(n+1,r2);
X	c.r+=rc.r*sdparams[10];
X	c.g+=rc.g*sdparams[11];
X	c.b+=rc.b*sdparams[12];
X	return c;
X}
X
Xstruct color glass(n,r,i)
Xint n;
Xstruct ray r;
Xstruct intersect i;
X{int j;
Xstruct color c,rc;
Xdouble l,x,y,z,sp;
Xstruct ray r2;
Xstruct vector e,li;
Xdouble *sdparams;
X	sdparams=objects[i.obj].sdparams;
X	c.r=sdparams[0];
X	c.g=sdparams[1];
X	c.b=sdparams[2];
X	x=r.a.x+i.t*r.l.x;
X	y=r.a.y+i.t*r.l.y;
X	z=r.a.z+i.t*r.l.z;
X	e.x=x-objects[0].center.x;
X	e.y=y-objects[0].center.y;
X	e.z=z-objects[0].center.z;
X	l=1/sqrt(e.x*e.x+e.y*e.y+e.z*e.z);
X	e.x*=l;
X	e.y*=l;
X	e.z*=l;
X	for(j=1;j<lights+1;j++)
X	{	r2.a=objects[j].center;
X		li.x=x-r2.a.x;
X		li.y=y-r2.a.y;
X		li.z=z-r2.a.z;
X		l=1/sqrt(li.x*li.x+li.y*li.y+li.z*li.z);
X		li.x*=l;
X		li.y*=l;
X		li.z*=l;
X		r2.obj=j;
X		r2.l=li;
X		l= -(i.n.x*li.x+i.n.y*li.y+i.n.z*li.z);
X		if(l>0&&!shadow(i.obj,r2))
X		{	c.r+=l*sdparams[3];
X			c.g+=l*sdparams[4];
X			c.b+=l*sdparams[5];
X			sp=phong(i.n,e,li,sdparams[9]);
X			c.r+=sp*sdparams[6];
X			c.g+=sp*sdparams[7];
X			c.b+=sp*sdparams[8];
X		}
X	}
X	r2.obj=i.obj;
X	r2.a.x=x;
X	r2.a.y=y;
X	r2.a.z=z;
X	r2.l=reflect(i.n,r.l);
X	rc=trace(n+1,r2);
X	c.r+=rc.r*sdparams[10];
X	c.g+=rc.g*sdparams[11];
X	c.b+=rc.b*sdparams[12];
X	r2.l=refract(i.n,r.l,sdparams[16]);
X	if(r2.l.x==0&&r2.l.y==0&&r2.l.z==0)
X		return c;
X	rc=trace(n+1,r2);
X	c.r+=rc.r*sdparams[13];
X	c.g+=rc.g*sdparams[14];
X	c.b+=rc.b*sdparams[15];
X	return c;
X}
EOF26639
sed 's/^X//' <<'EOF26640' >shape.c
X#include <math.h>
X
X#include "ray.h"
X
Xvoid initsphere();
Xstruct intersect sphere();
X
Xinitshape()
X{	nshapetab=1;
X	shapetab=nalloc(struct shapetab,nshapetab);
X	shapetab[0].shapeinitfunc=initsphere;
X	shapetab[0].shapefunc=sphere;
X	shapetab[0].nparams=4;
X}
X
Xvoid initsphere(o)
Xstruct object *o;
X{	o->center.x=o->spparams[0];
X	o->center.y=o->spparams[1];
X	o->center.z=o->spparams[2];
X	o->radius=o->spparams[3];
X}
X
Xstruct intersect sphere(r,obj)
Xstruct ray r;
Xregister int obj;
X{struct vector v,n;
Xregister double *p,b,c,d,t,t0,t1;
Xstruct intersect i;
X	i.obj=0;
X	p=objects[obj].spparams;
X	v.x=r.a.x-p[0];
X	v.y=r.a.y-p[1];
X	v.z=r.a.z-p[2];
X	b=r.l.x*v.x+r.l.y*v.y+r.l.z*v.z;
X	c=v.x*v.x+v.y*v.y+v.z*v.z-p[3]*p[3];
X	d=b*b-c;
X	if(d<0)
X		return i;
X	d=sqrt(d);
X	if(b>=0)
X	{	t0= -b-d;
X		t1=c/(-b-d);
X	}
X	else
X	{	t0=c/(-b+d);
X		t1= -b+d;
X	}
X	if(t0>MINT)
X		t=t0;
X	else if(t1>MINT)
X		t=t1;
X	else
X		return i;
X	i.obj=obj;
X	i.t=t;
X	n.x=r.a.x+r.l.x*t-p[0];
X	n.y=r.a.y+r.l.y*t-p[1];
X	n.z=r.a.z+r.l.z*t-p[2];
X	t=sqrt(n.x*n.x+n.y*n.y+n.z*n.z);
X	i.n.x=n.x/t;
X	i.n.y=n.y/t;
X	i.n.z=n.z/t;
X	return i;
X}
EOF26640
sed 's/^X//' <<'EOF26641' >trace.c
X#include "ray.h"
X
Xstruct color black={0,0,0},white={1.0,1.0,1.0};
X
Xstruct color trace(n,r)
Xregister int n;
Xstruct ray r;
X{struct intersect i;
X	if(n>=MAXLEVEL)
X		return bgcolor(r);
X	raycount++;
X	i=intersect(r);
X	if(i.obj>0)
X	{	if(sflag)
X			return (*objects[i.obj].shade)(n,r,i);
X		else
X			return white;
X	}
X	else
X	{	if(sflag)
X			return bgcolor(r);
X		else
X			return black;
X	}
X}
EOF26641
sed 's/^X//' <<'EOF26642' >test.v
Xf 109.5
Xl	10 10 -5
Xl	-20 0 -5
Xo 0 1
X	-5 -5 10 3
X	0.2 0 0		0.2 0 0		0.6 0.0 0.0 0.01
Xo 0 1
X	-3 5 10 2
X	0 0.2 0		0 0.2 0		0.0 0.6 0.0 0.01
Xo 0 1
X	6 -2 10 2
X	0 0 0.2		0 0 0.2		0.0 0.0 0.6 0.01
Xo 0 1
X	0 -4 20 7
X	0.2 0.2 0	0.2 0.2 0	0.6 0.6 0.6 0.01
Xo 0 2
X	5 5 15 5
X	0 0 0		0 0 0		0.9 0.9 0.9 0.01	0.9 0.9 0.9
Xo 0 3
X	-1 1 7 3
X	0 0 0		0 0 0		0.2 0.2 0.2 0.01	0.1 0.1 0.1
X	0.9 0.9 0.9 1.2
Xe
EOF26642
sed 's/^X//' <<'EOF26643' >random.c
X#include <stdio.h>
X#include <math.h>
X
Xdouble rnd()
X{	return (random()&0xffff)/65536.0;
X}
X
Xmain(ac,av)
Xint ac;
Xchar **av;
X{int i,j,n;
Xdouble x,y,z,r,cr,d,l;
X	if(ac!=3&&ac!=4)
X	{	fprintf(stderr,"Usage: %s number radious [seed]\n",
X			av[0]);
X		exit(1);
X	}
X	n=atoi(av[1]);
X	r=atof(av[2]);
X	if(ac==4)
X		srandom(atoi(av[3]));
X	if(n<0||r<=0)
X	{	fprintf(stderr,"Usage: %s number radious [seed]\n",
X			av[0]);
X		exit(1);
X	}
X	for(i=0;i<n;i++)
X	{	while(1)
X		{	x=rnd()*2-1;
X			y=rnd()*2-1;
X			z=rnd()*2-1;
X			l=x*x+y*y+z*z;
X			if(l<=(1-r)*(1-r))
X				break;;
X		}
X		printf("o 0 0\n");
X		printf("\t%g %g %g %g\n",x,y,z+2,r);
X		printf("\t0.2 0.2 0.2\t0.4 0.4 0.4\n");
X	}
X	exit(0);
X}
EOF26643
exit