Mercurial > cgi-bin > hgwebdir.cgi > PR > Applications > SSR > SSR__C-Ray__Bench
diff c-ray-mt.c @ 0:11a4bcadac2a
Initial SSR version
| author | Merten Sach <msach@mailbox.tu-berlin.de> |
|---|---|
| date | Thu, 22 Sep 2011 14:16:25 +0200 |
| parents | |
| children | b6c9e5f46e98 |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/c-ray-mt.c Thu Sep 22 14:16:25 2011 +0200 1.3 @@ -0,0 +1,704 @@ 1.4 +/* c-ray-mt - a simple multithreaded raytracing filter. 1.5 + * Copyright (C) 2006 John Tsiombikas <nuclear@siggraph.org> 1.6 + * 1.7 + * You are free to use, modify and redistribute this program under the 1.8 + * terms of the GNU General Public License v2 or (at your option) later. 1.9 + * see "http://www.gnu.org/licenses/gpl.txt" for details. 1.10 + * --------------------------------------------------------------------- 1.11 + * Usage: 1.12 + * compile: just type make 1.13 + * (add any arch-specific optimizations for your compiler in CFLAGS first) 1.14 + * run: cat scene | ./c-ray-mt [-t num-threads] >foo.ppm 1.15 + * (on broken systems such as windows try: c-ray-mt -i scene -o foo.ppm) 1.16 + * enjoy: display foo.ppm 1.17 + * (with imagemagick, or use your favorite image viewer) 1.18 + * --------------------------------------------------------------------- 1.19 + * Scene file format: 1.20 + * # sphere (many) 1.21 + * s x y z rad r g b shininess reflectivity 1.22 + * # light (many) 1.23 + * l x y z 1.24 + * # camera (one) 1.25 + * c x y z fov tx ty tz 1.26 + * --------------------------------------------------------------------- 1.27 + */ 1.28 +#include <stdio.h> 1.29 +#include <stdlib.h> 1.30 +#include <string.h> 1.31 +#include <math.h> 1.32 +#include <ctype.h> 1.33 +#include <errno.h> 1.34 +#include <pthread.h> 1.35 +#include "SSR_lib/SSR.h" 1.36 + 1.37 +#define VER_MAJOR 1 1.38 +#define VER_MINOR 1 1.39 +#define VER_STR "c-ray-mt v%d.%d\n" 1.40 + 1.41 +#if !defined(unix) && !defined(__unix__) 1.42 +#ifdef __MACH__ 1.43 +#define unix 1 1.44 +#define __unix__ 1 1.45 +#endif /* __MACH__ */ 1.46 +#endif /* unix */ 1.47 + 1.48 +/* find the appropriate way to define explicitly sized types */ 1.49 +/* for C99 or GNU libc (also mach's libc) we can use stdint.h */ 1.50 +#if (__STDC_VERSION__ >= 199900) || defined(__GLIBC__) || defined(__MACH__) 1.51 +#include <stdint.h> 1.52 +#elif defined(unix) || defined(__unix__) /* some UNIX systems have them in sys/types.h */ 1.53 +#include <sys/types.h> 1.54 +#elif defined(__WIN32__) || defined(WIN32) /* the nameless one */ 1.55 +typedef unsigned __int8 uint8_t; 1.56 +typedef unsigned __int32 uint32_t; 1.57 +#endif /* sized type detection */ 1.58 + 1.59 +struct vec3 { 1.60 + double x, y, z; 1.61 +}; 1.62 + 1.63 +struct ray { 1.64 + struct vec3 orig, dir; 1.65 +}; 1.66 + 1.67 +struct material { 1.68 + struct vec3 col; /* color */ 1.69 + double spow; /* specular power */ 1.70 + double refl; /* reflection intensity */ 1.71 +}; 1.72 + 1.73 +struct sphere { 1.74 + struct vec3 pos; 1.75 + double rad; 1.76 + struct material mat; 1.77 + struct sphere *next; 1.78 +}; 1.79 + 1.80 +struct spoint { 1.81 + struct vec3 pos, normal, vref; /* position, normal and view reflection */ 1.82 + double dist; /* parametric distance of intersection along the ray */ 1.83 +}; 1.84 + 1.85 +struct camera { 1.86 + struct vec3 pos, targ; 1.87 + double fov; 1.88 +}; 1.89 + 1.90 +struct procr_data { 1.91 + VirtProcr *VP; 1.92 + VirtProcr *parentVP; 1.93 + int sl_start, sl_count; 1.94 + uint32_t *pixels; 1.95 +}; 1.96 +typedef struct procr_data procr_data; 1.97 + 1.98 +void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples); 1.99 +struct vec3 trace(struct ray ray, int depth); 1.100 +struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth); 1.101 +struct vec3 reflect(struct vec3 v, struct vec3 n); 1.102 +struct vec3 cross_product(struct vec3 v1, struct vec3 v2); 1.103 +struct ray get_primary_ray(int x, int y, int sample); 1.104 +struct vec3 get_sample_pos(int x, int y, int sample); 1.105 +struct vec3 jitter(int x, int y, int s); 1.106 +int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp); 1.107 +void load_scene(FILE *fp); 1.108 +unsigned long get_msec(void); 1.109 + 1.110 +void thread_func(void *tdata, VirtProcr *VProc); 1.111 + 1.112 +#define MAX_LIGHTS 16 /* maximum number of lights */ 1.113 +#define RAY_MAG 1000.0 /* trace rays of this magnitude */ 1.114 +#define MAX_RAY_DEPTH 5 /* raytrace recursion limit */ 1.115 +#define FOV 0.78539816 /* field of view in rads (pi/4) */ 1.116 +#define HALF_FOV (FOV * 0.5) 1.117 +#define ERR_MARGIN 1e-6 /* an arbitrary error margin to avoid surface acne */ 1.118 + 1.119 +/* bit-shift ammount for packing each color into a 32bit uint */ 1.120 +#ifdef LITTLE_ENDIAN 1.121 +#define RSHIFT 16 1.122 +#define BSHIFT 0 1.123 +#else /* big endian */ 1.124 +#define RSHIFT 0 1.125 +#define BSHIFT 16 1.126 +#endif /* endianess */ 1.127 +#define GSHIFT 8 /* this is the same in both byte orders */ 1.128 + 1.129 +/* some helpful macros... */ 1.130 +#define SQ(x) ((x) * (x)) 1.131 +#define MAX(a, b) ((a) > (b) ? (a) : (b)) 1.132 +#define MIN(a, b) ((a) < (b) ? (a) : (b)) 1.133 +#define DOT(a, b) ((a).x * (b).x + (a).y * (b).y + (a).z * (b).z) 1.134 +#define NORMALIZE(a) do {\ 1.135 + double len = sqrt(DOT(a, a));\ 1.136 + (a).x /= len; (a).y /= len; (a).z /= len;\ 1.137 +} while(0); 1.138 + 1.139 +//SSR Message Types 1.140 +#define WORK_START 1 1.141 +#define WORK_END 2 1.142 + 1.143 +/* global state */ 1.144 +int xres = 800; 1.145 +int yres = 600; 1.146 +int rays_per_pixel = 1; 1.147 +double aspect = 1.333333; 1.148 +struct sphere *obj_list; 1.149 +struct vec3 lights[MAX_LIGHTS]; 1.150 +int lnum = 0; 1.151 +struct camera cam; 1.152 + 1.153 +int thread_num = 1; 1.154 +struct procr_data *procrs; 1.155 + 1.156 +volatile int end = 0; 1.157 +volatile int start = 0; 1.158 +int32 end_mutex, end_cond; 1.159 +int32 start_cond, start_mutex; 1.160 + 1.161 +#define NRAN 1024 1.162 +#define MASK (NRAN - 1) 1.163 +struct vec3 urand[NRAN]; 1.164 +int irand[NRAN]; 1.165 + 1.166 +unsigned long rend_time, start_time; 1.167 + 1.168 +const char *usage = { 1.169 + "Usage: c-ray-mt [options]\n" 1.170 + " Reads a scene file from stdin, writes the image to stdout, and stats to stderr.\n\n" 1.171 + "Options:\n" 1.172 + " -t <num> how many threads to use (default: 1)\n" 1.173 + " -s WxH where W is the width and H the height of the image\n" 1.174 + " -r <rays> shoot <rays> rays per pixel (antialiasing)\n" 1.175 + " -i <file> read from <file> instead of stdin\n" 1.176 + " -o <file> write to <file> instead of stdout\n" 1.177 + " -h this help screen\n\n" 1.178 +}; 1.179 + 1.180 +char __ProgrammName[] = "c-ray"; 1.181 +char __DataSet[255]; 1.182 + 1.183 + 1.184 +void raytrace(void *pixels, VirtProcr *Vprocr); 1.185 + 1.186 +int main(int argc, char **argv) { 1.187 + int i; 1.188 + uint32_t *pixels; 1.189 + FILE *infile = stdin, *outfile = stdout; 1.190 + 1.191 + for(i=1; i<argc; i++) { 1.192 + if(argv[i][0] == '-' && argv[i][2] == 0) { 1.193 + char *sep; 1.194 + switch(argv[i][1]) { 1.195 + case 't': 1.196 + if(!isdigit(argv[++i][0])) { 1.197 + fprintf(stderr, "-t mus be followed by the number of worker threads to spawn\n"); 1.198 + return EXIT_FAILURE; 1.199 + } 1.200 + thread_num = atoi(argv[i]); 1.201 + if(!thread_num) { 1.202 + fprintf(stderr, "invalid number of threads specified: %d\n", thread_num); 1.203 + return EXIT_FAILURE; 1.204 + } 1.205 + break; 1.206 + 1.207 + case 's': 1.208 + if(!isdigit(argv[++i][0]) || !(sep = strchr(argv[i], 'x')) || !isdigit(*(sep + 1))) { 1.209 + fputs("-s must be followed by something like \"640x480\"\n", stderr); 1.210 + return EXIT_FAILURE; 1.211 + } 1.212 + xres = atoi(argv[i]); 1.213 + yres = atoi(sep + 1); 1.214 + aspect = (double)xres / (double)yres; 1.215 + break; 1.216 + 1.217 + case 'i': 1.218 + if(!(infile = fopen(argv[++i], "rb"))) { 1.219 + fprintf(stderr, "failed to open input file %s: %s\n", argv[i], strerror(errno)); 1.220 + return EXIT_FAILURE; 1.221 + } 1.222 + break; 1.223 + 1.224 + case 'o': 1.225 + if(!(outfile = fopen(argv[++i], "wb"))) { 1.226 + fprintf(stderr, "failed to open output file %s: %s\n", argv[i], strerror(errno)); 1.227 + return EXIT_FAILURE; 1.228 + } 1.229 + break; 1.230 + 1.231 + case 'r': 1.232 + if(!isdigit(argv[++i][0])) { 1.233 + fputs("-r must be followed by a number (rays per pixel)\n", stderr); 1.234 + return EXIT_FAILURE; 1.235 + } 1.236 + rays_per_pixel = atoi(argv[i]); 1.237 + break; 1.238 + 1.239 + case 'h': 1.240 + fputs(usage, stdout); 1.241 + return 0; 1.242 + 1.243 + default: 1.244 + fprintf(stderr, "unrecognized argument: %s\n", argv[i]); 1.245 + fputs(usage, stderr); 1.246 + return EXIT_FAILURE; 1.247 + } 1.248 + } else { 1.249 + fprintf(stderr, "unrecognized argument: %s\n", argv[i]); 1.250 + fputs(usage, stderr); 1.251 + return EXIT_FAILURE; 1.252 + } 1.253 + } 1.254 + 1.255 + 1.256 + if(!(pixels = malloc(xres * yres * sizeof *pixels))) { 1.257 + perror("pixel buffer allocation failed"); 1.258 + return EXIT_FAILURE; 1.259 + } 1.260 + load_scene(infile); 1.261 + 1.262 + //This is the transition to the VMS runtime 1.263 + SSR__create_seed_procr_and_do_work(raytrace, (void*)pixels); 1.264 + 1.265 + /* output statistics to stderr */ 1.266 + fprintf(stderr, "Rendering took: %lu seconds (%lu milliseconds)\n", rend_time / 1000, rend_time); 1.267 + 1.268 + /* output the image */ 1.269 + fprintf(outfile, "P6\n%d %d\n255\n", xres, yres); 1.270 + for(i=0; i<xres * yres; i++) { 1.271 + fputc((pixels[i] >> RSHIFT) & 0xff, outfile); 1.272 + fputc((pixels[i] >> GSHIFT) & 0xff, outfile); 1.273 + fputc((pixels[i] >> BSHIFT) & 0xff, outfile); 1.274 + } 1.275 + fflush(outfile); 1.276 + 1.277 + if(infile != stdin) fclose(infile); 1.278 + if(outfile != stdout) fclose(outfile); 1.279 + 1.280 + struct sphere *walker = obj_list; 1.281 + while(walker) { 1.282 + struct sphere *tmp = walker; 1.283 + walker = walker->next; 1.284 + free(tmp); 1.285 + } 1.286 + free(pixels); 1.287 + return 0; 1.288 +} 1.289 + 1.290 +/* this is run after the VMS is set up*/ 1.291 +void raytrace(void *pixels, VirtProcr *VProc) 1.292 +{ 1.293 + int i; 1.294 + double sl, sl_per_procr; 1.295 + 1.296 + /* initialize the random number tables for the jitter */ 1.297 + for(i=0; i<NRAN; i++) urand[i].x = (double)rand() / RAND_MAX - 0.5; 1.298 + for(i=0; i<NRAN; i++) urand[i].y = (double)rand() / RAND_MAX - 0.5; 1.299 + for(i=0; i<NRAN; i++) irand[i] = (int)(NRAN * ((double)rand() / RAND_MAX)); 1.300 + 1.301 + if(thread_num > yres) { 1.302 + fprintf(stderr, "more threads than scanlines specified, reducing number of threads to %d\n", yres); 1.303 + thread_num = yres; 1.304 + } 1.305 + 1.306 + 1.307 + if(!(procrs = SSR__malloc_to(thread_num * sizeof(procr_data), VProc))) { 1.308 + perror("failed to allocate thread table"); 1.309 + exit(EXIT_FAILURE); 1.310 + } 1.311 + 1.312 + sl = 0.0; 1.313 + sl_per_procr = (double)yres / (double)thread_num; 1.314 + for(i=0; i<thread_num; i++) { 1.315 + procrs[i].sl_start = (int)sl; 1.316 + sl += sl_per_procr; 1.317 + procrs[i].sl_count = (int)sl - procrs[i].sl_start; 1.318 + procrs[i].pixels = (uint32_t*)pixels; 1.319 + procrs[i].parentVP = VProc; 1.320 + 1.321 + procrs[i].VP = SSR__create_procr_with((VirtProcrFnPtr)thread_func, 1.322 + (void*)(&procrs[i]), VProc); 1.323 + } 1.324 + 1.325 + procrs[thread_num - 1].sl_count = yres - procrs[thread_num - 1].sl_start; 1.326 + 1.327 + fprintf(stderr, VER_STR, VER_MAJOR, VER_MINOR); 1.328 + 1.329 + // start worker threads 1.330 + //printf("start of worker thread (%d)\n", VProc->procrID); 1.331 + start_time = get_msec(); 1.332 + for(i=0; i<thread_num; i++) 1.333 + SSR__send_of_type_to(VProc, NULL, WORK_START, procrs[i].VP); 1.334 + 1.335 + //printf("wait for worker (%d)\n", VProc->procrID); 1.336 + for(i=0; i<thread_num; i++) 1.337 + SSR__receive_type_to(WORK_END, VProc); 1.338 + 1.339 + rend_time = get_msec() - start_time; 1.340 + 1.341 + SSR__free(procrs,VProc); 1.342 + SSR__dissipate_procr(VProc); 1.343 +} 1.344 + 1.345 +/* render a frame of xsz/ysz dimensions into the provided framebuffer */ 1.346 +void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples) { 1.347 + int i, s; 1.348 + double rcp_samples = 1.0 / (double)samples; 1.349 + 1.350 + for(i=0; i<xsz; i++) { 1.351 + double r, g, b; 1.352 + r = g = b = 0.0; 1.353 + 1.354 + for(s=0; s<samples; s++) { 1.355 + struct vec3 col = trace(get_primary_ray(i, sl, s), 0); 1.356 + r += col.x; 1.357 + g += col.y; 1.358 + b += col.z; 1.359 + } 1.360 + 1.361 + r = r * rcp_samples; 1.362 + g = g * rcp_samples; 1.363 + b = b * rcp_samples; 1.364 + 1.365 + fb[sl * xsz + i] = ((uint32_t)(MIN(r, 1.0) * 255.0) & 0xff) << RSHIFT | 1.366 + ((uint32_t)(MIN(g, 1.0) * 255.0) & 0xff) << GSHIFT | 1.367 + ((uint32_t)(MIN(b, 1.0) * 255.0) & 0xff) << BSHIFT; 1.368 + } 1.369 +} 1.370 + 1.371 +/* trace a ray throught the scene recursively (the recursion happens through 1.372 + * shade() to calculate reflection rays if necessary). 1.373 + */ 1.374 +struct vec3 trace(struct ray ray, int depth) { 1.375 + struct vec3 col; 1.376 + struct spoint sp, nearest_sp; 1.377 + struct sphere *nearest_obj = 0; 1.378 + struct sphere *iter = obj_list->next; 1.379 + 1.380 + /* if we reached the recursion limit, bail out */ 1.381 + if(depth >= MAX_RAY_DEPTH) { 1.382 + col.x = col.y = col.z = 0.0; 1.383 + return col; 1.384 + } 1.385 + 1.386 + /* find the nearest intersection ... */ 1.387 + while(iter) { 1.388 + if(ray_sphere(iter, ray, &sp)) { 1.389 + if(!nearest_obj || sp.dist < nearest_sp.dist) { 1.390 + nearest_obj = iter; 1.391 + nearest_sp = sp; 1.392 + } 1.393 + } 1.394 + iter = iter->next; 1.395 + } 1.396 + 1.397 + /* and perform shading calculations as needed by calling shade() */ 1.398 + if(nearest_obj) { 1.399 + col = shade(nearest_obj, &nearest_sp, depth); 1.400 + } else { 1.401 + col.x = col.y = col.z = 0.0; 1.402 + } 1.403 + 1.404 + return col; 1.405 +} 1.406 + 1.407 +/* Calculates direct illumination with the phong reflectance model. 1.408 + * Also handles reflections by calling trace again, if necessary. 1.409 + */ 1.410 +struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth) { 1.411 + int i; 1.412 + struct vec3 col = {0, 0, 0}; 1.413 + 1.414 + /* for all lights ... */ 1.415 + for(i=0; i<lnum; i++) { 1.416 + double ispec, idiff; 1.417 + struct vec3 ldir; 1.418 + struct ray shadow_ray; 1.419 + struct sphere *iter = obj_list->next; 1.420 + int in_shadow = 0; 1.421 + 1.422 + ldir.x = lights[i].x - sp->pos.x; 1.423 + ldir.y = lights[i].y - sp->pos.y; 1.424 + ldir.z = lights[i].z - sp->pos.z; 1.425 + 1.426 + shadow_ray.orig = sp->pos; 1.427 + shadow_ray.dir = ldir; 1.428 + 1.429 + /* shoot shadow rays to determine if we have a line of sight with the light */ 1.430 + while(iter) { 1.431 + if(ray_sphere(iter, shadow_ray, 0)) { 1.432 + in_shadow = 1; 1.433 + break; 1.434 + } 1.435 + iter = iter->next; 1.436 + } 1.437 + 1.438 + /* and if we're not in shadow, calculate direct illumination with the phong model. */ 1.439 + if(!in_shadow) { 1.440 + NORMALIZE(ldir); 1.441 + 1.442 + idiff = MAX(DOT(sp->normal, ldir), 0.0); 1.443 + ispec = obj->mat.spow > 0.0 ? pow(MAX(DOT(sp->vref, ldir), 0.0), obj->mat.spow) : 0.0; 1.444 + 1.445 + col.x += idiff * obj->mat.col.x + ispec; 1.446 + col.y += idiff * obj->mat.col.y + ispec; 1.447 + col.z += idiff * obj->mat.col.z + ispec; 1.448 + } 1.449 + } 1.450 + 1.451 + /* Also, if the object is reflective, spawn a reflection ray, and call trace() 1.452 + * to calculate the light arriving from the mirror direction. 1.453 + */ 1.454 + if(obj->mat.refl > 0.0) { 1.455 + struct ray ray; 1.456 + struct vec3 rcol; 1.457 + 1.458 + ray.orig = sp->pos; 1.459 + ray.dir = sp->vref; 1.460 + ray.dir.x *= RAY_MAG; 1.461 + ray.dir.y *= RAY_MAG; 1.462 + ray.dir.z *= RAY_MAG; 1.463 + 1.464 + rcol = trace(ray, depth + 1); 1.465 + col.x += rcol.x * obj->mat.refl; 1.466 + col.y += rcol.y * obj->mat.refl; 1.467 + col.z += rcol.z * obj->mat.refl; 1.468 + } 1.469 + 1.470 + return col; 1.471 +} 1.472 + 1.473 +/* calculate reflection vector */ 1.474 +struct vec3 reflect(struct vec3 v, struct vec3 n) { 1.475 + struct vec3 res; 1.476 + double dot = v.x * n.x + v.y * n.y + v.z * n.z; 1.477 + res.x = -(2.0 * dot * n.x - v.x); 1.478 + res.y = -(2.0 * dot * n.y - v.y); 1.479 + res.z = -(2.0 * dot * n.z - v.z); 1.480 + return res; 1.481 +} 1.482 + 1.483 +struct vec3 cross_product(struct vec3 v1, struct vec3 v2) { 1.484 + struct vec3 res; 1.485 + res.x = v1.y * v2.z - v1.z * v2.y; 1.486 + res.y = v1.z * v2.x - v1.x * v2.z; 1.487 + res.z = v1.x * v2.y - v1.y * v2.x; 1.488 + return res; 1.489 +} 1.490 + 1.491 +/* determine the primary ray corresponding to the specified pixel (x, y) */ 1.492 +struct ray get_primary_ray(int x, int y, int sample) { 1.493 + struct ray ray; 1.494 + float m[3][3]; 1.495 + struct vec3 i, j = {0, 1, 0}, k, dir, orig, foo; 1.496 + 1.497 + k.x = cam.targ.x - cam.pos.x; 1.498 + k.y = cam.targ.y - cam.pos.y; 1.499 + k.z = cam.targ.z - cam.pos.z; 1.500 + NORMALIZE(k); 1.501 + 1.502 + i = cross_product(j, k); 1.503 + j = cross_product(k, i); 1.504 + m[0][0] = i.x; m[0][1] = j.x; m[0][2] = k.x; 1.505 + m[1][0] = i.y; m[1][1] = j.y; m[1][2] = k.y; 1.506 + m[2][0] = i.z; m[2][1] = j.z; m[2][2] = k.z; 1.507 + 1.508 + ray.orig.x = ray.orig.y = ray.orig.z = 0.0; 1.509 + ray.dir = get_sample_pos(x, y, sample); 1.510 + ray.dir.z = 1.0 / HALF_FOV; 1.511 + ray.dir.x *= RAY_MAG; 1.512 + ray.dir.y *= RAY_MAG; 1.513 + ray.dir.z *= RAY_MAG; 1.514 + 1.515 + dir.x = ray.dir.x + ray.orig.x; 1.516 + dir.y = ray.dir.y + ray.orig.y; 1.517 + dir.z = ray.dir.z + ray.orig.z; 1.518 + foo.x = dir.x * m[0][0] + dir.y * m[0][1] + dir.z * m[0][2]; 1.519 + foo.y = dir.x * m[1][0] + dir.y * m[1][1] + dir.z * m[1][2]; 1.520 + foo.z = dir.x * m[2][0] + dir.y * m[2][1] + dir.z * m[2][2]; 1.521 + 1.522 + orig.x = ray.orig.x * m[0][0] + ray.orig.y * m[0][1] + ray.orig.z * m[0][2] + cam.pos.x; 1.523 + orig.y = ray.orig.x * m[1][0] + ray.orig.y * m[1][1] + ray.orig.z * m[1][2] + cam.pos.y; 1.524 + orig.z = ray.orig.x * m[2][0] + ray.orig.y * m[2][1] + ray.orig.z * m[2][2] + cam.pos.z; 1.525 + 1.526 + ray.orig = orig; 1.527 + ray.dir.x = foo.x + orig.x; 1.528 + ray.dir.y = foo.y + orig.y; 1.529 + ray.dir.z = foo.z + orig.z; 1.530 + 1.531 + return ray; 1.532 +} 1.533 + 1.534 + 1.535 +struct vec3 get_sample_pos(int x, int y, int sample) { 1.536 + struct vec3 pt; 1.537 + static double sf = 0.0; 1.538 + 1.539 + if(sf == 0.0) { 1.540 + sf = 1.5 / (double)xres; 1.541 + } 1.542 + 1.543 + pt.x = ((double)x / (double)xres) - 0.5; 1.544 + pt.y = -(((double)y / (double)yres) - 0.65) / aspect; 1.545 + 1.546 + if(sample) { 1.547 + struct vec3 jt = jitter(x, y, sample); 1.548 + pt.x += jt.x * sf; 1.549 + pt.y += jt.y * sf / aspect; 1.550 + } 1.551 + return pt; 1.552 +} 1.553 + 1.554 +/* jitter function taken from Graphics Gems I. */ 1.555 +struct vec3 jitter(int x, int y, int s) { 1.556 + struct vec3 pt; 1.557 + pt.x = urand[(x + (y << 2) + irand[(x + s) & MASK]) & MASK].x; 1.558 + pt.y = urand[(y + (x << 2) + irand[(y + s) & MASK]) & MASK].y; 1.559 + return pt; 1.560 +} 1.561 + 1.562 +/* Calculate ray-sphere intersection, and return {1, 0} to signify hit or no hit. 1.563 + * Also the surface point parameters like position, normal, etc are returned through 1.564 + * the sp pointer if it is not NULL. 1.565 + */ 1.566 +int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp) { 1.567 + double a, b, c, d, sqrt_d, t1, t2; 1.568 + 1.569 + a = SQ(ray.dir.x) + SQ(ray.dir.y) + SQ(ray.dir.z); 1.570 + b = 2.0 * ray.dir.x * (ray.orig.x - sph->pos.x) + 1.571 + 2.0 * ray.dir.y * (ray.orig.y - sph->pos.y) + 1.572 + 2.0 * ray.dir.z * (ray.orig.z - sph->pos.z); 1.573 + c = SQ(sph->pos.x) + SQ(sph->pos.y) + SQ(sph->pos.z) + 1.574 + SQ(ray.orig.x) + SQ(ray.orig.y) + SQ(ray.orig.z) + 1.575 + 2.0 * (-sph->pos.x * ray.orig.x - sph->pos.y * ray.orig.y - sph->pos.z * ray.orig.z) - SQ(sph->rad); 1.576 + 1.577 + if((d = SQ(b) - 4.0 * a * c) < 0.0) return 0; 1.578 + 1.579 + sqrt_d = sqrt(d); 1.580 + t1 = (-b + sqrt_d) / (2.0 * a); 1.581 + t2 = (-b - sqrt_d) / (2.0 * a); 1.582 + 1.583 + if((t1 < ERR_MARGIN && t2 < ERR_MARGIN) || (t1 > 1.0 && t2 > 1.0)) return 0; 1.584 + 1.585 + if(sp) { 1.586 + if(t1 < ERR_MARGIN) t1 = t2; 1.587 + if(t2 < ERR_MARGIN) t2 = t1; 1.588 + sp->dist = t1 < t2 ? t1 : t2; 1.589 + 1.590 + sp->pos.x = ray.orig.x + ray.dir.x * sp->dist; 1.591 + sp->pos.y = ray.orig.y + ray.dir.y * sp->dist; 1.592 + sp->pos.z = ray.orig.z + ray.dir.z * sp->dist; 1.593 + 1.594 + sp->normal.x = (sp->pos.x - sph->pos.x) / sph->rad; 1.595 + sp->normal.y = (sp->pos.y - sph->pos.y) / sph->rad; 1.596 + sp->normal.z = (sp->pos.z - sph->pos.z) / sph->rad; 1.597 + 1.598 + sp->vref = reflect(ray.dir, sp->normal); 1.599 + NORMALIZE(sp->vref); 1.600 + } 1.601 + return 1; 1.602 +} 1.603 + 1.604 +/* Load the scene from an extremely simple scene description file */ 1.605 +#define DELIM " \t\n" 1.606 +void load_scene(FILE *fp) { 1.607 + char line[256], *ptr, type; 1.608 + 1.609 + obj_list = malloc(sizeof(struct sphere)); 1.610 + obj_list->next = 0; 1.611 + 1.612 + while((ptr = fgets(line, 256, fp))) { 1.613 + int i; 1.614 + struct vec3 pos, col; 1.615 + double rad, spow, refl; 1.616 + 1.617 + while(*ptr == ' ' || *ptr == '\t') ptr++; 1.618 + if(*ptr == '#' || *ptr == '\n') continue; 1.619 + 1.620 + if(!(ptr = strtok(line, DELIM))) continue; 1.621 + type = *ptr; 1.622 + 1.623 + for(i=0; i<3; i++) { 1.624 + if(!(ptr = strtok(0, DELIM))) break; 1.625 + *((double*)&pos.x + i) = atof(ptr); 1.626 + } 1.627 + 1.628 + if(type == 'l') { 1.629 + lights[lnum++] = pos; 1.630 + continue; 1.631 + } 1.632 + 1.633 + if(!(ptr = strtok(0, DELIM))) continue; 1.634 + rad = atof(ptr); 1.635 + 1.636 + for(i=0; i<3; i++) { 1.637 + if(!(ptr = strtok(0, DELIM))) break; 1.638 + *((double*)&col.x + i) = atof(ptr); 1.639 + } 1.640 + 1.641 + if(type == 'c') { 1.642 + cam.pos = pos; 1.643 + cam.targ = col; 1.644 + cam.fov = rad; 1.645 + continue; 1.646 + } 1.647 + 1.648 + if(!(ptr = strtok(0, DELIM))) continue; 1.649 + spow = atof(ptr); 1.650 + 1.651 + if(!(ptr = strtok(0, DELIM))) continue; 1.652 + refl = atof(ptr); 1.653 + 1.654 + if(type == 's') { 1.655 + struct sphere *sph = malloc(sizeof *sph); 1.656 + sph->next = obj_list->next; 1.657 + obj_list->next = sph; 1.658 + 1.659 + sph->pos = pos; 1.660 + sph->rad = rad; 1.661 + sph->mat.col = col; 1.662 + sph->mat.spow = spow; 1.663 + sph->mat.refl = refl; 1.664 + } else { 1.665 + fprintf(stderr, "unknown type: %c\n", type); 1.666 + } 1.667 + } 1.668 +} 1.669 + 1.670 + 1.671 +/* provide a millisecond-resolution timer for each system */ 1.672 +#if defined(unix) || defined(__unix__) 1.673 +#include <time.h> 1.674 +#include <sys/time.h> 1.675 +unsigned long get_msec(void) { 1.676 + static struct timeval timeval, first_timeval; 1.677 + 1.678 + gettimeofday(&timeval, 0); 1.679 + if(first_timeval.tv_sec == 0) { 1.680 + first_timeval = timeval; 1.681 + return 0; 1.682 + } 1.683 + return (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000; 1.684 +} 1.685 +#elif defined(__WIN32__) || defined(WIN32) 1.686 +#include <windows.h> 1.687 +unsigned long get_msec(void) { 1.688 + return GetTickCount(); 1.689 +} 1.690 +#else 1.691 +#error "I don't know how to measure time on your platform" 1.692 +#endif 1.693 + 1.694 +void thread_func(void *tdata, VirtProcr *VProc) { 1.695 + int i; 1.696 + procr_data *td = (procr_data*)tdata; 1.697 + 1.698 + SSR__receive_type_to(WORK_START, VProc); 1.699 + 1.700 + for(i=0; i<td->sl_count; i++) { 1.701 + render_scanline(xres, yres, i + td->sl_start, td->pixels, rays_per_pixel); 1.702 + } 1.703 + 1.704 + SSR__send_of_type_to(VProc, NULL, WORK_END, td->parentVP); 1.705 + 1.706 + SSR__dissipate_procr(VProc); 1.707 +}