msach@0: /* c-ray-mt - a simple multithreaded raytracing filter.
msach@0:  * Copyright (C) 2006 John Tsiombikas <nuclear@siggraph.org>
msach@0:  *
msach@0:  * You are free to use, modify and redistribute this program under the
msach@0:  * terms of the GNU General Public License v2 or (at your option) later.
msach@0:  * see "http://www.gnu.org/licenses/gpl.txt" for details.
msach@0:  * ---------------------------------------------------------------------
msach@0:  * Usage:
msach@0:  *   compile:  just type make
msach@0:  *              (add any arch-specific optimizations for your compiler in CFLAGS first)
msach@0:  *       run:  cat scene | ./c-ray-mt [-t num-threads] >foo.ppm
msach@0:  *              (on broken systems such as windows try: c-ray-mt -i scene -o foo.ppm)
msach@0:  *     enjoy:  display foo.ppm
msach@0:  *              (with imagemagick, or use your favorite image viewer)
msach@0:  * ---------------------------------------------------------------------
msach@0:  * Scene file format:
msach@0:  *   # sphere (many)
msach@0:  *   s  x y z  rad   r g b   shininess   reflectivity
msach@0:  *   # light (many)
msach@0:  *   l  x y z
msach@0:  *   # camera (one)
msach@0:  *   c  x y z  fov   tx ty tz
msach@0:  * ---------------------------------------------------------------------
msach@0:  */
msach@0: #include <stdio.h>
msach@0: #include <stdlib.h>
msach@0: #include <string.h>
msach@0: #include <math.h>
msach@0: #include <ctype.h>
msach@0: #include <errno.h>
msach@0: #include <pthread.h>
msach@0: #include "VPThread_lib/VPThread.h"
msach@0: 
msach@0: #define VER_MAJOR	1
msach@0: #define VER_MINOR	1
msach@0: #define VER_STR		"c-ray-mt v%d.%d\n"
msach@0: 
msach@0: #if !defined(unix) && !defined(__unix__)
msach@0: #ifdef __MACH__
msach@0: #define unix		1
msach@0: #define __unix__	1
msach@0: #endif	/* __MACH__ */
msach@0: #endif	/* unix */
msach@0: 
msach@0: /* find the appropriate way to define explicitly sized types */
msach@0: /* for C99 or GNU libc (also mach's libc) we can use stdint.h */
msach@0: #if (__STDC_VERSION__ >= 199900) || defined(__GLIBC__) || defined(__MACH__)
msach@0: #include <stdint.h>
msach@0: #elif defined(unix) || defined(__unix__)	/* some UNIX systems have them in sys/types.h */
msach@0: #include <sys/types.h>
msach@0: #elif defined(__WIN32__) || defined(WIN32)	/* the nameless one */
msach@0: typedef unsigned __int8 uint8_t;
msach@0: typedef unsigned __int32 uint32_t;
msach@0: #endif	/* sized type detection */
msach@0: 
msach@0: struct vec3 {
msach@0: 	double x, y, z;
msach@0: };
msach@0: 
msach@0: struct ray {
msach@0: 	struct vec3 orig, dir;
msach@0: };
msach@0: 
msach@0: struct material {
msach@0: 	struct vec3 col;	/* color */
msach@0: 	double spow;		/* specular power */
msach@0: 	double refl;		/* reflection intensity */
msach@0: };
msach@0: 
msach@0: struct sphere {
msach@0: 	struct vec3 pos;
msach@0: 	double rad;
msach@0: 	struct material mat;
msach@0: 	struct sphere *next;
msach@0: };
msach@0: 
msach@0: struct spoint {
msach@0: 	struct vec3 pos, normal, vref;	/* position, normal and view reflection */
msach@0: 	double dist;		/* parametric distance of intersection along the ray */
msach@0: };
msach@0: 
msach@0: struct camera {
msach@0: 	struct vec3 pos, targ;
msach@0: 	double fov;
msach@0: };
msach@0: 
msach@0: struct thread_data {
msach@1: 	VirtProcr *VP;
msach@0: 	int sl_start, sl_count;
msach@0: 	uint32_t *pixels;
msach@0: };
msach@1: typedef struct thread_data thread_data;
msach@0: 
msach@0: void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples);
msach@0: struct vec3 trace(struct ray ray, int depth);
msach@0: struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth);
msach@0: struct vec3 reflect(struct vec3 v, struct vec3 n);
msach@0: struct vec3 cross_product(struct vec3 v1, struct vec3 v2);
msach@0: struct ray get_primary_ray(int x, int y, int sample);
msach@0: struct vec3 get_sample_pos(int x, int y, int sample);
msach@0: struct vec3 jitter(int x, int y, int s);
msach@0: int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp);
msach@0: void load_scene(FILE *fp);
msach@0: unsigned long get_msec(void);
msach@0: 
msach@1: void thread_func(void *tdata, VirtProcr *VProc);
msach@0: 
msach@0: #define MAX_LIGHTS		16				/* maximum number of lights */
msach@0: #define RAY_MAG			1000.0			/* trace rays of this magnitude */
msach@0: #define MAX_RAY_DEPTH	5				/* raytrace recursion limit */
msach@0: #define FOV				0.78539816		/* field of view in rads (pi/4) */
msach@0: #define HALF_FOV		(FOV * 0.5)
msach@0: #define ERR_MARGIN		1e-6			/* an arbitrary error margin to avoid surface acne */
msach@0: 
msach@0: /* bit-shift ammount for packing each color into a 32bit uint */
msach@0: #ifdef LITTLE_ENDIAN
msach@0: #define RSHIFT	16
msach@0: #define BSHIFT	0
msach@0: #else	/* big endian */
msach@0: #define RSHIFT	0
msach@0: #define BSHIFT	16
msach@0: #endif	/* endianess */
msach@0: #define GSHIFT	8	/* this is the same in both byte orders */
msach@0: 
msach@0: /* some helpful macros... */
msach@0: #define SQ(x)		((x) * (x))
msach@0: #define MAX(a, b)	((a) > (b) ? (a) : (b))
msach@0: #define MIN(a, b)	((a) < (b) ? (a) : (b))
msach@0: #define DOT(a, b)	((a).x * (b).x + (a).y * (b).y + (a).z * (b).z)
msach@0: #define NORMALIZE(a)  do {\
msach@0: 	double len = sqrt(DOT(a, a));\
msach@0: 	(a).x /= len; (a).y /= len; (a).z /= len;\
msach@0: } while(0);
msach@0: 
msach@0: /* global state */
msach@0: int xres = 800;
msach@0: int yres = 600;
msach@0: int rays_per_pixel = 1;
msach@0: double aspect = 1.333333;
msach@0: struct sphere *obj_list;
msach@0: struct vec3 lights[MAX_LIGHTS];
msach@0: int lnum = 0;
msach@0: struct camera cam;
msach@0: 
msach@0: int thread_num = 1;
msach@0: struct thread_data *threads;
msach@0: 
msach@1: volatile int end = 0;
msach@1: volatile int start = 0;
msach@1: int32 end_mutex, end_cond;
msach@1: int32 start_cond, start_mutex;
msach@0: 
msach@0: #define NRAN	1024
msach@0: #define MASK	(NRAN - 1)
msach@0: struct vec3 urand[NRAN];
msach@0: int irand[NRAN];
msach@0: 
msach@0: unsigned long rend_time, start_time;
msach@0: 
msach@0: const char *usage = {
msach@0: 	"Usage: c-ray-mt [options]\n"
msach@0: 	"  Reads a scene file from stdin, writes the image to stdout, and stats to stderr.\n\n"
msach@0: 	"Options:\n"
msach@0: 	"  -t <num>   how many threads to use (default: 1)\n"
msach@0: 	"  -s WxH     where W is the width and H the height of the image\n"
msach@0: 	"  -r <rays>  shoot <rays> rays per pixel (antialiasing)\n"
msach@0: 	"  -i <file>  read from <file> instead of stdin\n"
msach@0: 	"  -o <file>  write to <file> instead of stdout\n"
msach@0: 	"  -h         this help screen\n\n"
msach@0: };
msach@0: 
msach@1: char __ProgrammName[] = "c-ray";
msach@1: char __DataSet[255];
msach@1: 
msach@1: 
msach@1: void raytrace(void *pixels, VirtProcr *Vprocr);
msach@0: 
msach@0: int main(int argc, char **argv) {
msach@0: 	int i;
msach@0: 	uint32_t *pixels;
msach@0: 	FILE *infile = stdin, *outfile = stdout;
msach@0: 
msach@0: 	for(i=1; i<argc; i++) {
msach@0: 		if(argv[i][0] == '-' && argv[i][2] == 0) {
msach@0: 			char *sep;
msach@0: 			switch(argv[i][1]) {
msach@0: 			case 't':
msach@0: 				if(!isdigit(argv[++i][0])) {
msach@0: 					fprintf(stderr, "-t mus be followed by the number of worker threads to spawn\n");
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				thread_num = atoi(argv[i]);
msach@0: 				if(!thread_num) {
msach@0: 					fprintf(stderr, "invalid number of threads specified: %d\n", thread_num);
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				break;
msach@0: 					
msach@0: 			case 's':
msach@0: 				if(!isdigit(argv[++i][0]) || !(sep = strchr(argv[i], 'x')) || !isdigit(*(sep + 1))) {
msach@0: 					fputs("-s must be followed by something like \"640x480\"\n", stderr);
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				xres = atoi(argv[i]);
msach@0: 				yres = atoi(sep + 1);
msach@0: 				aspect = (double)xres / (double)yres;
msach@0: 				break;
msach@0: 
msach@0: 			case 'i':
msach@0: 				if(!(infile = fopen(argv[++i], "rb"))) {
msach@0: 					fprintf(stderr, "failed to open input file %s: %s\n", argv[i], strerror(errno));
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				break;
msach@0: 
msach@0: 			case 'o':
msach@0: 				if(!(outfile = fopen(argv[++i], "wb"))) {
msach@0: 					fprintf(stderr, "failed to open output file %s: %s\n", argv[i], strerror(errno));
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				break;
msach@0: 
msach@0: 			case 'r':
msach@0: 				if(!isdigit(argv[++i][0])) {
msach@0: 					fputs("-r must be followed by a number (rays per pixel)\n", stderr);
msach@0: 					return EXIT_FAILURE;
msach@0: 				}
msach@0: 				rays_per_pixel = atoi(argv[i]);
msach@0: 				break;
msach@0: 
msach@0: 			case 'h':
msach@0: 				fputs(usage, stdout);
msach@0: 				return 0;
msach@0: 				
msach@0: 			default:
msach@0: 				fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
msach@0: 				fputs(usage, stderr);
msach@0: 				return EXIT_FAILURE;
msach@0: 			}
msach@0: 		} else {
msach@0: 			fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
msach@0: 			fputs(usage, stderr);
msach@0: 			return EXIT_FAILURE;
msach@0: 		}
msach@0: 	}
msach@0: 
msach@0:         
msach@0:         if(!(pixels = malloc(xres * yres * sizeof *pixels))) {
msach@0: 		perror("pixel buffer allocation failed");
msach@0: 		return EXIT_FAILURE;
msach@0: 	}
msach@0: 	load_scene(infile);
msach@0:         
msach@1:         //This is the transition to the VMS runtime
msach@1:         VPThread__create_seed_procr_and_do_work(raytrace, (void*)pixels);
msach@0: 	
msach@0: 	/* output statistics to stderr */
msach@0: 	fprintf(stderr, "Rendering took: %lu seconds (%lu milliseconds)\n", rend_time / 1000, rend_time);
msach@0: 
msach@0: 	/* output the image */
msach@0: 	fprintf(outfile, "P6\n%d %d\n255\n", xres, yres);
msach@0: 	for(i=0; i<xres * yres; i++) {
msach@0: 		fputc((pixels[i] >> RSHIFT) & 0xff, outfile);
msach@0: 		fputc((pixels[i] >> GSHIFT) & 0xff, outfile);
msach@0: 		fputc((pixels[i] >> BSHIFT) & 0xff, outfile);
msach@0: 	}
msach@0: 	fflush(outfile);
msach@0: 
msach@0: 	if(infile != stdin) fclose(infile);
msach@0: 	if(outfile != stdout) fclose(outfile);
msach@0: 
msach@0: 	struct sphere *walker = obj_list;
msach@0: 	while(walker) {
msach@0: 		struct sphere *tmp = walker;
msach@0: 		walker = walker->next;
msach@0: 		free(tmp);
msach@0: 	}
msach@0: 	free(pixels);
msach@0: 	return 0;
msach@0: }
msach@0: 
msach@0: /* this is run after the VMS is set up*/
msach@1: void raytrace(void *pixels, VirtProcr *VProc)
msach@0: {
msach@0:     int i;
msach@0:     double sl, sl_per_thread;
msach@0:     
msach@0:     /* initialize the random number tables for the jitter */
msach@0:     for(i=0; i<NRAN; i++) urand[i].x = (double)rand() / RAND_MAX - 0.5;
msach@0:     for(i=0; i<NRAN; i++) urand[i].y = (double)rand() / RAND_MAX - 0.5;
msach@0:     for(i=0; i<NRAN; i++) irand[i] = (int)(NRAN * ((double)rand() / RAND_MAX));
msach@0: 
msach@0:     if(thread_num > yres) {
msach@0:             fprintf(stderr, "more threads than scanlines specified, reducing number of threads to %d\n", yres);
msach@0:             thread_num = yres;
msach@0:     }
msach@0: 
msach@1:     
msach@1:     if(!(threads = VPThread__malloc(thread_num * sizeof(thread_data), VProc))) {
msach@0:             perror("failed to allocate thread table");
msach@0:             exit(EXIT_FAILURE);
msach@0:     }
msach@1:     
msach@1:     end_mutex = VPThread__make_mutex(VProc);
msach@1:     end_cond  = VPThread__make_cond(end_mutex, VProc);
msach@1:     start_mutex = VPThread__make_mutex(VProc);
msach@1:     start_cond  = VPThread__make_cond(start_mutex, VProc);    
msach@1:     
msach@0:     sl = 0.0;
msach@0:     sl_per_thread = (double)yres / (double)thread_num;
msach@0:     for(i=0; i<thread_num; i++) {
msach@0:             threads[i].sl_start = (int)sl;
msach@0:             sl += sl_per_thread;
msach@0:             threads[i].sl_count = (int)sl - threads[i].sl_start;
msach@1:             threads[i].pixels = (uint32_t*)pixels;
msach@0: 
msach@1:             threads[i].VP = 
msach@1:                     VPThread__create_thread((VirtProcrFnPtr)thread_func,
msach@1:                                    (void*)(&threads[i]), VProc);  
msach@0:     }
msach@1:     
msach@0:     threads[thread_num - 1].sl_count = yres - threads[thread_num - 1].sl_start;
msach@1:     
msach@0:     fprintf(stderr, VER_STR, VER_MAJOR, VER_MINOR);
msach@1:     
msach@1:     // start worker threads
msach@1:     //printf("start of worker thread (%d)\n", VProc->procrID);      
msach@1:     VPThread__mutex_lock(start_mutex, VProc);
msach@0:     start_time = get_msec();
msach@0:     start = 1;
msach@1:     for(i=0; i<thread_num; i++)
msach@1:         VPThread__cond_signal(start_cond, VProc);
msach@1:     VPThread__mutex_unlock(start_mutex, VProc);
msach@1:     
msach@1:     //printf("wait for worker (%d)\n", VProc->procrID);      
msach@1:     VPThread__mutex_lock(end_mutex, VProc);
msach@1:     while(end < thread_num)
msach@1:         VPThread__cond_wait(end_cond, VProc);
msach@1:     VPThread__mutex_unlock(end_mutex, VProc);
msach@1:     
msach@0:     rend_time = get_msec() - start_time;
msach@1:     
msach@1:     VPThread__free(threads,VProc);
msach@1:     VPThread__dissipate_thread(VProc);
msach@0: }
msach@0: 
msach@0: /* render a frame of xsz/ysz dimensions into the provided framebuffer */
msach@0: void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples) {
msach@0: 	int i, s;
msach@0: 	double rcp_samples = 1.0 / (double)samples;
msach@0: 
msach@0: 	for(i=0; i<xsz; i++) {
msach@0: 		double r, g, b;
msach@0: 		r = g = b = 0.0;
msach@0: 			
msach@0: 		for(s=0; s<samples; s++) {
msach@0: 			struct vec3 col = trace(get_primary_ray(i, sl, s), 0);
msach@0: 			r += col.x;
msach@0: 			g += col.y;
msach@0: 			b += col.z;
msach@0: 		}
msach@0: 
msach@0: 		r = r * rcp_samples;
msach@0: 		g = g * rcp_samples;
msach@0: 		b = b * rcp_samples;
msach@0: 			
msach@0: 		fb[sl * xsz + i] = ((uint32_t)(MIN(r, 1.0) * 255.0) & 0xff) << RSHIFT |
msach@0: 							((uint32_t)(MIN(g, 1.0) * 255.0) & 0xff) << GSHIFT |
msach@0: 							((uint32_t)(MIN(b, 1.0) * 255.0) & 0xff) << BSHIFT;
msach@0: 	}
msach@0: }
msach@0: 
msach@0: /* trace a ray throught the scene recursively (the recursion happens through
msach@0:  * shade() to calculate reflection rays if necessary).
msach@0:  */
msach@0: struct vec3 trace(struct ray ray, int depth) {
msach@0: 	struct vec3 col;
msach@0: 	struct spoint sp, nearest_sp;
msach@0: 	struct sphere *nearest_obj = 0;
msach@0: 	struct sphere *iter = obj_list->next;
msach@0: 
msach@0: 	/* if we reached the recursion limit, bail out */
msach@0: 	if(depth >= MAX_RAY_DEPTH) {
msach@0: 		col.x = col.y = col.z = 0.0;
msach@0: 		return col;
msach@0: 	}
msach@0: 	
msach@0: 	/* find the nearest intersection ... */
msach@0: 	while(iter) {
msach@0: 		if(ray_sphere(iter, ray, &sp)) {
msach@0: 			if(!nearest_obj || sp.dist < nearest_sp.dist) {
msach@0: 				nearest_obj = iter;
msach@0: 				nearest_sp = sp;
msach@0: 			}
msach@0: 		}
msach@0: 		iter = iter->next;
msach@0: 	}
msach@0: 
msach@0: 	/* and perform shading calculations as needed by calling shade() */
msach@0: 	if(nearest_obj) {
msach@0: 		col = shade(nearest_obj, &nearest_sp, depth);
msach@0: 	} else {
msach@0: 		col.x = col.y = col.z = 0.0;
msach@0: 	}
msach@0: 
msach@0: 	return col;
msach@0: }
msach@0: 
msach@0: /* Calculates direct illumination with the phong reflectance model.
msach@0:  * Also handles reflections by calling trace again, if necessary.
msach@0:  */
msach@0: struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth) {
msach@0: 	int i;
msach@0: 	struct vec3 col = {0, 0, 0};
msach@0: 
msach@0: 	/* for all lights ... */
msach@0: 	for(i=0; i<lnum; i++) {
msach@0: 		double ispec, idiff;
msach@0: 		struct vec3 ldir;
msach@0: 		struct ray shadow_ray;
msach@0: 		struct sphere *iter = obj_list->next;
msach@0: 		int in_shadow = 0;
msach@0: 
msach@0: 		ldir.x = lights[i].x - sp->pos.x;
msach@0: 		ldir.y = lights[i].y - sp->pos.y;
msach@0: 		ldir.z = lights[i].z - sp->pos.z;
msach@0: 
msach@0: 		shadow_ray.orig = sp->pos;
msach@0: 		shadow_ray.dir = ldir;
msach@0: 
msach@0: 		/* shoot shadow rays to determine if we have a line of sight with the light */
msach@0: 		while(iter) {
msach@0: 			if(ray_sphere(iter, shadow_ray, 0)) {
msach@0: 				in_shadow = 1;
msach@0: 				break;
msach@0: 			}
msach@0: 			iter = iter->next;
msach@0: 		}
msach@0: 
msach@0: 		/* and if we're not in shadow, calculate direct illumination with the phong model. */
msach@0: 		if(!in_shadow) {
msach@0: 			NORMALIZE(ldir);
msach@0: 
msach@0: 			idiff = MAX(DOT(sp->normal, ldir), 0.0);
msach@0: 			ispec = obj->mat.spow > 0.0 ? pow(MAX(DOT(sp->vref, ldir), 0.0), obj->mat.spow) : 0.0;
msach@0: 
msach@0: 			col.x += idiff * obj->mat.col.x + ispec;
msach@0: 			col.y += idiff * obj->mat.col.y + ispec;
msach@0: 			col.z += idiff * obj->mat.col.z + ispec;
msach@0: 		}
msach@0: 	}
msach@0: 
msach@0: 	/* Also, if the object is reflective, spawn a reflection ray, and call trace()
msach@0: 	 * to calculate the light arriving from the mirror direction.
msach@0: 	 */
msach@0: 	if(obj->mat.refl > 0.0) {
msach@0: 		struct ray ray;
msach@0: 		struct vec3 rcol;
msach@0: 
msach@0: 		ray.orig = sp->pos;
msach@0: 		ray.dir = sp->vref;
msach@0: 		ray.dir.x *= RAY_MAG;
msach@0: 		ray.dir.y *= RAY_MAG;
msach@0: 		ray.dir.z *= RAY_MAG;
msach@0: 
msach@0: 		rcol = trace(ray, depth + 1);
msach@0: 		col.x += rcol.x * obj->mat.refl;
msach@0: 		col.y += rcol.y * obj->mat.refl;
msach@0: 		col.z += rcol.z * obj->mat.refl;
msach@0: 	}
msach@0: 
msach@0: 	return col;
msach@0: }
msach@0: 
msach@0: /* calculate reflection vector */
msach@0: struct vec3 reflect(struct vec3 v, struct vec3 n) {
msach@0: 	struct vec3 res;
msach@0: 	double dot = v.x * n.x + v.y * n.y + v.z * n.z;
msach@0: 	res.x = -(2.0 * dot * n.x - v.x);
msach@0: 	res.y = -(2.0 * dot * n.y - v.y);
msach@0: 	res.z = -(2.0 * dot * n.z - v.z);
msach@0: 	return res;
msach@0: }
msach@0: 
msach@0: struct vec3 cross_product(struct vec3 v1, struct vec3 v2) {
msach@0: 	struct vec3 res;
msach@0: 	res.x = v1.y * v2.z - v1.z * v2.y;
msach@0: 	res.y = v1.z * v2.x - v1.x * v2.z;
msach@0: 	res.z = v1.x * v2.y - v1.y * v2.x;
msach@0: 	return res;
msach@0: }
msach@0: 
msach@0: /* determine the primary ray corresponding to the specified pixel (x, y) */
msach@0: struct ray get_primary_ray(int x, int y, int sample) {
msach@0: 	struct ray ray;
msach@0: 	float m[3][3];
msach@0: 	struct vec3 i, j = {0, 1, 0}, k, dir, orig, foo;
msach@0: 
msach@0: 	k.x = cam.targ.x - cam.pos.x;
msach@0: 	k.y = cam.targ.y - cam.pos.y;
msach@0: 	k.z = cam.targ.z - cam.pos.z;
msach@0: 	NORMALIZE(k);
msach@0: 
msach@0: 	i = cross_product(j, k);
msach@0: 	j = cross_product(k, i);
msach@0: 	m[0][0] = i.x; m[0][1] = j.x; m[0][2] = k.x;
msach@0: 	m[1][0] = i.y; m[1][1] = j.y; m[1][2] = k.y;
msach@0: 	m[2][0] = i.z; m[2][1] = j.z; m[2][2] = k.z;
msach@0: 	
msach@0: 	ray.orig.x = ray.orig.y = ray.orig.z = 0.0;
msach@0: 	ray.dir = get_sample_pos(x, y, sample);
msach@0: 	ray.dir.z = 1.0 / HALF_FOV;
msach@0: 	ray.dir.x *= RAY_MAG;
msach@0: 	ray.dir.y *= RAY_MAG;
msach@0: 	ray.dir.z *= RAY_MAG;
msach@0: 	
msach@0: 	dir.x = ray.dir.x + ray.orig.x;
msach@0: 	dir.y = ray.dir.y + ray.orig.y;
msach@0: 	dir.z = ray.dir.z + ray.orig.z;
msach@0: 	foo.x = dir.x * m[0][0] + dir.y * m[0][1] + dir.z * m[0][2];
msach@0: 	foo.y = dir.x * m[1][0] + dir.y * m[1][1] + dir.z * m[1][2];
msach@0: 	foo.z = dir.x * m[2][0] + dir.y * m[2][1] + dir.z * m[2][2];
msach@0: 
msach@0: 	orig.x = ray.orig.x * m[0][0] + ray.orig.y * m[0][1] + ray.orig.z * m[0][2] + cam.pos.x;
msach@0: 	orig.y = ray.orig.x * m[1][0] + ray.orig.y * m[1][1] + ray.orig.z * m[1][2] + cam.pos.y;
msach@0: 	orig.z = ray.orig.x * m[2][0] + ray.orig.y * m[2][1] + ray.orig.z * m[2][2] + cam.pos.z;
msach@0: 
msach@0: 	ray.orig = orig;
msach@0: 	ray.dir.x = foo.x + orig.x;
msach@0: 	ray.dir.y = foo.y + orig.y;
msach@0: 	ray.dir.z = foo.z + orig.z;
msach@0: 	
msach@0: 	return ray;
msach@0: }
msach@0: 
msach@0: 
msach@0: struct vec3 get_sample_pos(int x, int y, int sample) {
msach@0: 	struct vec3 pt;
msach@0: 	static double sf = 0.0;
msach@0: 
msach@0: 	if(sf == 0.0) {
msach@0: 		sf = 1.5 / (double)xres;
msach@0: 	}
msach@0: 
msach@0: 	pt.x = ((double)x / (double)xres) - 0.5;
msach@0: 	pt.y = -(((double)y / (double)yres) - 0.65) / aspect;
msach@0: 
msach@0: 	if(sample) {
msach@0: 		struct vec3 jt = jitter(x, y, sample);
msach@0: 		pt.x += jt.x * sf;
msach@0: 		pt.y += jt.y * sf / aspect;
msach@0: 	}
msach@0: 	return pt;
msach@0: }
msach@0: 
msach@0: /* jitter function taken from Graphics Gems I. */
msach@0: struct vec3 jitter(int x, int y, int s) {
msach@0: 	struct vec3 pt;
msach@0: 	pt.x = urand[(x + (y << 2) + irand[(x + s) & MASK]) & MASK].x;
msach@0: 	pt.y = urand[(y + (x << 2) + irand[(y + s) & MASK]) & MASK].y;
msach@0: 	return pt;
msach@0: }
msach@0: 
msach@0: /* Calculate ray-sphere intersection, and return {1, 0} to signify hit or no hit.
msach@0:  * Also the surface point parameters like position, normal, etc are returned through
msach@0:  * the sp pointer if it is not NULL.
msach@0:  */
msach@0: int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp) {
msach@0: 	double a, b, c, d, sqrt_d, t1, t2;
msach@0: 	
msach@0: 	a = SQ(ray.dir.x) + SQ(ray.dir.y) + SQ(ray.dir.z);
msach@0: 	b = 2.0 * ray.dir.x * (ray.orig.x - sph->pos.x) +
msach@0: 				2.0 * ray.dir.y * (ray.orig.y - sph->pos.y) +
msach@0: 				2.0 * ray.dir.z * (ray.orig.z - sph->pos.z);
msach@0: 	c = SQ(sph->pos.x) + SQ(sph->pos.y) + SQ(sph->pos.z) +
msach@0: 				SQ(ray.orig.x) + SQ(ray.orig.y) + SQ(ray.orig.z) +
msach@0: 				2.0 * (-sph->pos.x * ray.orig.x - sph->pos.y * ray.orig.y - sph->pos.z * ray.orig.z) - SQ(sph->rad);
msach@0: 	
msach@0: 	if((d = SQ(b) - 4.0 * a * c) < 0.0) return 0;
msach@0: 
msach@0: 	sqrt_d = sqrt(d);
msach@0: 	t1 = (-b + sqrt_d) / (2.0 * a);
msach@0: 	t2 = (-b - sqrt_d) / (2.0 * a);
msach@0: 
msach@0: 	if((t1 < ERR_MARGIN && t2 < ERR_MARGIN) || (t1 > 1.0 && t2 > 1.0)) return 0;
msach@0: 
msach@0: 	if(sp) {
msach@0: 		if(t1 < ERR_MARGIN) t1 = t2;
msach@0: 		if(t2 < ERR_MARGIN) t2 = t1;
msach@0: 		sp->dist = t1 < t2 ? t1 : t2;
msach@0: 		
msach@0: 		sp->pos.x = ray.orig.x + ray.dir.x * sp->dist;
msach@0: 		sp->pos.y = ray.orig.y + ray.dir.y * sp->dist;
msach@0: 		sp->pos.z = ray.orig.z + ray.dir.z * sp->dist;
msach@0: 		
msach@0: 		sp->normal.x = (sp->pos.x - sph->pos.x) / sph->rad;
msach@0: 		sp->normal.y = (sp->pos.y - sph->pos.y) / sph->rad;
msach@0: 		sp->normal.z = (sp->pos.z - sph->pos.z) / sph->rad;
msach@0: 
msach@0: 		sp->vref = reflect(ray.dir, sp->normal);
msach@0: 		NORMALIZE(sp->vref);
msach@0: 	}
msach@0: 	return 1;
msach@0: }
msach@0: 
msach@0: /* Load the scene from an extremely simple scene description file */
msach@0: #define DELIM	" \t\n"
msach@0: void load_scene(FILE *fp) {
msach@0: 	char line[256], *ptr, type;
msach@0: 
msach@0: 	obj_list = malloc(sizeof(struct sphere));
msach@0: 	obj_list->next = 0;
msach@0: 	
msach@0: 	while((ptr = fgets(line, 256, fp))) {
msach@0: 		int i;
msach@0: 		struct vec3 pos, col;
msach@0: 		double rad, spow, refl;
msach@0: 		
msach@0: 		while(*ptr == ' ' || *ptr == '\t') ptr++;
msach@0: 		if(*ptr == '#' || *ptr == '\n') continue;
msach@0: 
msach@0: 		if(!(ptr = strtok(line, DELIM))) continue;
msach@0: 		type = *ptr;
msach@0: 		
msach@0: 		for(i=0; i<3; i++) {
msach@0: 			if(!(ptr = strtok(0, DELIM))) break;
msach@0: 			*((double*)&pos.x + i) = atof(ptr);
msach@0: 		}
msach@0: 
msach@0: 		if(type == 'l') {
msach@0: 			lights[lnum++] = pos;
msach@0: 			continue;
msach@0: 		}
msach@0: 
msach@0: 		if(!(ptr = strtok(0, DELIM))) continue;
msach@0: 		rad = atof(ptr);
msach@0: 
msach@0: 		for(i=0; i<3; i++) {
msach@0: 			if(!(ptr = strtok(0, DELIM))) break;
msach@0: 			*((double*)&col.x + i) = atof(ptr);
msach@0: 		}
msach@0: 
msach@0: 		if(type == 'c') {
msach@0: 			cam.pos = pos;
msach@0: 			cam.targ = col;
msach@0: 			cam.fov = rad;
msach@0: 			continue;
msach@0: 		}
msach@0: 
msach@0: 		if(!(ptr = strtok(0, DELIM))) continue;
msach@0: 		spow = atof(ptr);
msach@0: 
msach@0: 		if(!(ptr = strtok(0, DELIM))) continue;
msach@0: 		refl = atof(ptr);
msach@0: 
msach@0: 		if(type == 's') {
msach@0: 			struct sphere *sph = malloc(sizeof *sph);
msach@0: 			sph->next = obj_list->next;
msach@0: 			obj_list->next = sph;
msach@0: 
msach@0: 			sph->pos = pos;
msach@0: 			sph->rad = rad;
msach@0: 			sph->mat.col = col;
msach@0: 			sph->mat.spow = spow;
msach@0: 			sph->mat.refl = refl;
msach@0: 		} else {
msach@0: 			fprintf(stderr, "unknown type: %c\n", type);
msach@0: 		}
msach@0: 	}
msach@0: }
msach@0: 
msach@0: 
msach@0: /* provide a millisecond-resolution timer for each system */
msach@0: #if defined(unix) || defined(__unix__)
msach@0: #include <time.h>
msach@0: #include <sys/time.h>
msach@0: unsigned long get_msec(void) {
msach@0: 	static struct timeval timeval, first_timeval;
msach@0: 	
msach@0: 	gettimeofday(&timeval, 0);
msach@0: 	if(first_timeval.tv_sec == 0) {
msach@0: 		first_timeval = timeval;
msach@0: 		return 0;
msach@0: 	}
msach@0: 	return (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000;
msach@0: }
msach@0: #elif defined(__WIN32__) || defined(WIN32)
msach@0: #include <windows.h>
msach@0: unsigned long get_msec(void) {
msach@0: 	return GetTickCount();
msach@0: }
msach@0: #else
msach@0: #error "I don't know how to measure time on your platform"
msach@0: #endif
msach@0: 
msach@1: void thread_func(void *tdata, VirtProcr *VProc) {
msach@0: 	int i;
msach@0: 	struct thread_data *td = (struct thread_data*)tdata;
msach@0: 
msach@1:         VPThread__mutex_lock(start_mutex, VProc);
msach@1:         while(!start)
msach@1:             VPThread__cond_wait(start_cond, VProc);
msach@1:         VPThread__mutex_unlock(start_mutex, VProc);        
msach@1:         
msach@0: 	for(i=0; i<td->sl_count; i++) {
msach@0: 		render_scanline(xres, yres, i + td->sl_start, td->pixels, rays_per_pixel);
msach@0: 	}
msach@1:         
msach@1:         VPThread__mutex_lock(end_mutex, VProc);
msach@1:         end++;
msach@1:         VPThread__cond_signal(end_cond, VProc);
msach@1:         VPThread__mutex_unlock(end_mutex, VProc);   
msach@0: 
msach@1: 	VPThread__dissipate_thread(VProc);
msach@0: }