Mercurial > cgi-bin > hgwebdir.cgi > PR > Applications > SSR > SSR__C-Ray__Bench

annotate c-ray-mt.c @ 2:8e3accb31ff3

find changesets by author, revision, files, or words in the commit message
typo corrected
author Merten Sach <msach@mailbox.tu-berlin.de>
date Wed, 28 Sep 2011 16:14:54 +0200
parents b6c9e5f46e98
children
rev   line source
msach@0 1 /* c-ray-mt - a simple multithreaded raytracing filter.
msach@0 2 * Copyright (C) 2006 John Tsiombikas <nuclear@siggraph.org>
msach@0 3 *
msach@0 4 * You are free to use, modify and redistribute this program under the
msach@0 5 * terms of the GNU General Public License v2 or (at your option) later.
msach@0 6 * see "http://www.gnu.org/licenses/gpl.txt" for details.
msach@0 7 * ---------------------------------------------------------------------
msach@0 8 * Usage:
msach@0 9 * compile: just type make
msach@0 10 * (add any arch-specific optimizations for your compiler in CFLAGS first)
msach@0 11 * run: cat scene | ./c-ray-mt [-t num-threads] >foo.ppm
msach@0 12 * (on broken systems such as windows try: c-ray-mt -i scene -o foo.ppm)
msach@0 13 * enjoy: display foo.ppm
msach@0 14 * (with imagemagick, or use your favorite image viewer)
msach@0 15 * ---------------------------------------------------------------------
msach@0 16 * Scene file format:
msach@0 17 * # sphere (many)
msach@0 18 * s x y z rad r g b shininess reflectivity
msach@0 19 * # light (many)
msach@0 20 * l x y z
msach@0 21 * # camera (one)
msach@0 22 * c x y z fov tx ty tz
msach@0 23 * ---------------------------------------------------------------------
msach@0 24 */
msach@0 25 #include <stdio.h>
msach@0 26 #include <stdlib.h>
msach@0 27 #include <string.h>
msach@0 28 #include <math.h>
msach@0 29 #include <ctype.h>
msach@0 30 #include <errno.h>
msach@0 31 #include <pthread.h>
msach@0 32 #include "SSR_lib/SSR.h"
msach@0 33
msach@0 34 #define VER_MAJOR 1
msach@0 35 #define VER_MINOR 1
msach@0 36 #define VER_STR "c-ray-mt v%d.%d\n"
msach@0 37
msach@0 38 #if !defined(unix) && !defined(__unix__)
msach@0 39 #ifdef __MACH__
msach@0 40 #define unix 1
msach@0 41 #define __unix__ 1
msach@0 42 #endif /* __MACH__ */
msach@0 43 #endif /* unix */
msach@0 44
msach@0 45 /* find the appropriate way to define explicitly sized types */
msach@0 46 /* for C99 or GNU libc (also mach's libc) we can use stdint.h */
msach@0 47 #if (__STDC_VERSION__ >= 199900) || defined(__GLIBC__) || defined(__MACH__)
msach@0 48 #include <stdint.h>
msach@0 49 #elif defined(unix) || defined(__unix__) /* some UNIX systems have them in sys/types.h */
msach@0 50 #include <sys/types.h>
msach@0 51 #elif defined(__WIN32__) || defined(WIN32) /* the nameless one */
msach@0 52 typedef unsigned __int8 uint8_t;
msach@0 53 typedef unsigned __int32 uint32_t;
msach@0 54 #endif /* sized type detection */
msach@0 55
msach@0 56 struct vec3 {
msach@0 57 double x, y, z;
msach@0 58 };
msach@0 59
msach@0 60 struct ray {
msach@0 61 struct vec3 orig, dir;
msach@0 62 };
msach@0 63
msach@0 64 struct material {
msach@0 65 struct vec3 col; /* color */
msach@0 66 double spow; /* specular power */
msach@0 67 double refl; /* reflection intensity */
msach@0 68 };
msach@0 69
msach@0 70 struct sphere {
msach@0 71 struct vec3 pos;
msach@0 72 double rad;
msach@0 73 struct material mat;
msach@0 74 struct sphere *next;
msach@0 75 };
msach@0 76
msach@0 77 struct spoint {
msach@0 78 struct vec3 pos, normal, vref; /* position, normal and view reflection */
msach@0 79 double dist; /* parametric distance of intersection along the ray */
msach@0 80 };
msach@0 81
msach@0 82 struct camera {
msach@0 83 struct vec3 pos, targ;
msach@0 84 double fov;
msach@0 85 };
msach@0 86
msach@0 87 struct procr_data {
msach@0 88 VirtProcr *VP;
msach@0 89 VirtProcr *parentVP;
msach@0 90 int sl_start, sl_count;
msach@0 91 uint32_t *pixels;
msach@0 92 };
msach@0 93 typedef struct procr_data procr_data;
msach@0 94
msach@0 95 void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples);
msach@0 96 struct vec3 trace(struct ray ray, int depth);
msach@0 97 struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth);
msach@0 98 struct vec3 reflect(struct vec3 v, struct vec3 n);
msach@0 99 struct vec3 cross_product(struct vec3 v1, struct vec3 v2);
msach@0 100 struct ray get_primary_ray(int x, int y, int sample);
msach@0 101 struct vec3 get_sample_pos(int x, int y, int sample);
msach@0 102 struct vec3 jitter(int x, int y, int s);
msach@0 103 int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp);
msach@0 104 void load_scene(FILE *fp);
msach@0 105 unsigned long get_msec(void);
msach@0 106
msach@0 107 void thread_func(void *tdata, VirtProcr *VProc);
msach@0 108
msach@0 109 #define MAX_LIGHTS 16 /* maximum number of lights */
msach@0 110 #define RAY_MAG 1000.0 /* trace rays of this magnitude */
msach@0 111 #define MAX_RAY_DEPTH 5 /* raytrace recursion limit */
msach@0 112 #define FOV 0.78539816 /* field of view in rads (pi/4) */
msach@0 113 #define HALF_FOV (FOV * 0.5)
msach@0 114 #define ERR_MARGIN 1e-6 /* an arbitrary error margin to avoid surface acne */
msach@0 115
msach@0 116 /* bit-shift ammount for packing each color into a 32bit uint */
msach@0 117 #ifdef LITTLE_ENDIAN
msach@0 118 #define RSHIFT 16
msach@0 119 #define BSHIFT 0
msach@0 120 #else /* big endian */
msach@0 121 #define RSHIFT 0
msach@0 122 #define BSHIFT 16
msach@0 123 #endif /* endianess */
msach@0 124 #define GSHIFT 8 /* this is the same in both byte orders */
msach@0 125
msach@0 126 /* some helpful macros... */
msach@0 127 #define SQ(x) ((x) * (x))
msach@0 128 #define MAX(a, b) ((a) > (b) ? (a) : (b))
msach@0 129 #define MIN(a, b) ((a) < (b) ? (a) : (b))
msach@0 130 #define DOT(a, b) ((a).x * (b).x + (a).y * (b).y + (a).z * (b).z)
msach@0 131 #define NORMALIZE(a) do {\
msach@0 132 double len = sqrt(DOT(a, a));\
msach@0 133 (a).x /= len; (a).y /= len; (a).z /= len;\
msach@0 134 } while(0);
msach@0 135
msach@0 136 //SSR Message Types
msach@0 137 #define WORK_START 1
msach@0 138 #define WORK_END 2
msach@0 139
msach@0 140 /* global state */
msach@0 141 int xres = 800;
msach@0 142 int yres = 600;
msach@0 143 int rays_per_pixel = 1;
msach@0 144 double aspect = 1.333333;
msach@0 145 struct sphere *obj_list;
msach@0 146 struct vec3 lights[MAX_LIGHTS];
msach@0 147 int lnum = 0;
msach@0 148 struct camera cam;
msach@0 149
msach@0 150 int thread_num = 1;
msach@0 151 struct procr_data *procrs;
msach@0 152
msach@0 153 volatile int end = 0;
msach@0 154 volatile int start = 0;
msach@0 155 int32 end_mutex, end_cond;
msach@0 156 int32 start_cond, start_mutex;
msach@0 157
msach@0 158 #define NRAN 1024
msach@0 159 #define MASK (NRAN - 1)
msach@0 160 struct vec3 urand[NRAN];
msach@0 161 int irand[NRAN];
msach@0 162
msach@0 163 unsigned long rend_time, start_time;
msach@0 164
msach@0 165 const char *usage = {
msach@0 166 "Usage: c-ray-mt [options]\n"
msach@0 167 " Reads a scene file from stdin, writes the image to stdout, and stats to stderr.\n\n"
msach@0 168 "Options:\n"
msach@0 169 " -t <num> how many threads to use (default: 1)\n"
msach@0 170 " -s WxH where W is the width and H the height of the image\n"
msach@0 171 " -r <rays> shoot <rays> rays per pixel (antialiasing)\n"
msach@0 172 " -i <file> read from <file> instead of stdin\n"
msach@0 173 " -o <file> write to <file> instead of stdout\n"
msach@0 174 " -h this help screen\n\n"
msach@0 175 };
msach@0 176
msach@0 177 char __ProgrammName[] = "c-ray";
msach@0 178 char __DataSet[255];
msach@0 179
msach@0 180
msach@0 181 void raytrace(void *pixels, VirtProcr *Vprocr);
msach@0 182
msach@0 183 int main(int argc, char **argv) {
msach@0 184 int i;
msach@0 185 uint32_t *pixels;
msach@0 186 FILE *infile = stdin, *outfile = stdout;
msach@0 187
msach@0 188 for(i=1; i<argc; i++) {
msach@0 189 if(argv[i][0] == '-' && argv[i][2] == 0) {
msach@0 190 char *sep;
msach@0 191 switch(argv[i][1]) {
msach@0 192 case 't':
msach@0 193 if(!isdigit(argv[++i][0])) {
msach@0 194 fprintf(stderr, "-t mus be followed by the number of worker threads to spawn\n");
msach@0 195 return EXIT_FAILURE;
msach@0 196 }
msach@0 197 thread_num = atoi(argv[i]);
msach@0 198 if(!thread_num) {
msach@0 199 fprintf(stderr, "invalid number of threads specified: %d\n", thread_num);
msach@0 200 return EXIT_FAILURE;
msach@0 201 }
msach@0 202 break;
msach@0 203
msach@0 204 case 's':
msach@0 205 if(!isdigit(argv[++i][0]) || !(sep = strchr(argv[i], 'x')) || !isdigit(*(sep + 1))) {
msach@0 206 fputs("-s must be followed by something like \"640x480\"\n", stderr);
msach@0 207 return EXIT_FAILURE;
msach@0 208 }
msach@0 209 xres = atoi(argv[i]);
msach@0 210 yres = atoi(sep + 1);
msach@0 211 aspect = (double)xres / (double)yres;
msach@0 212 break;
msach@0 213
msach@0 214 case 'i':
msach@0 215 if(!(infile = fopen(argv[++i], "rb"))) {
msach@0 216 fprintf(stderr, "failed to open input file %s: %s\n", argv[i], strerror(errno));
msach@0 217 return EXIT_FAILURE;
msach@0 218 }
msach@0 219 break;
msach@0 220
msach@0 221 case 'o':
msach@0 222 if(!(outfile = fopen(argv[++i], "wb"))) {
msach@0 223 fprintf(stderr, "failed to open output file %s: %s\n", argv[i], strerror(errno));
msach@0 224 return EXIT_FAILURE;
msach@0 225 }
msach@0 226 break;
msach@0 227
msach@0 228 case 'r':
msach@0 229 if(!isdigit(argv[++i][0])) {
msach@0 230 fputs("-r must be followed by a number (rays per pixel)\n", stderr);
msach@0 231 return EXIT_FAILURE;
msach@0 232 }
msach@0 233 rays_per_pixel = atoi(argv[i]);
msach@0 234 break;
msach@0 235
msach@0 236 case 'h':
msach@0 237 fputs(usage, stdout);
msach@0 238 return 0;
msach@0 239
msach@0 240 default:
msach@0 241 fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
msach@0 242 fputs(usage, stderr);
msach@0 243 return EXIT_FAILURE;
msach@0 244 }
msach@0 245 } else {
msach@0 246 fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
msach@0 247 fputs(usage, stderr);
msach@0 248 return EXIT_FAILURE;
msach@0 249 }
msach@0 250 }
msach@1 251
msach@2 252 snprintf(__DataSet,255,"file: %s\nsize: %dx%d\nrays per pixel: %d\nthreads: %d\n",
msach@1 253 infile, xres, yres, rays_per_pixel, thread_num);
msach@0 254
msach@0 255 if(!(pixels = malloc(xres * yres * sizeof *pixels))) {
msach@0 256 perror("pixel buffer allocation failed");
msach@0 257 return EXIT_FAILURE;
msach@0 258 }
msach@0 259 load_scene(infile);
msach@0 260
msach@0 261 //This is the transition to the VMS runtime
msach@0 262 SSR__create_seed_procr_and_do_work(raytrace, (void*)pixels);
msach@0 263
msach@0 264 /* output statistics to stderr */
msach@0 265 fprintf(stderr, "Rendering took: %lu seconds (%lu milliseconds)\n", rend_time / 1000, rend_time);
msach@0 266
msach@0 267 /* output the image */
msach@0 268 fprintf(outfile, "P6\n%d %d\n255\n", xres, yres);
msach@0 269 for(i=0; i<xres * yres; i++) {
msach@0 270 fputc((pixels[i] >> RSHIFT) & 0xff, outfile);
msach@0 271 fputc((pixels[i] >> GSHIFT) & 0xff, outfile);
msach@0 272 fputc((pixels[i] >> BSHIFT) & 0xff, outfile);
msach@0 273 }
msach@0 274 fflush(outfile);
msach@0 275
msach@0 276 if(infile != stdin) fclose(infile);
msach@0 277 if(outfile != stdout) fclose(outfile);
msach@0 278
msach@0 279 struct sphere *walker = obj_list;
msach@0 280 while(walker) {
msach@0 281 struct sphere *tmp = walker;
msach@0 282 walker = walker->next;
msach@0 283 free(tmp);
msach@0 284 }
msach@0 285 free(pixels);
msach@0 286 return 0;
msach@0 287 }
msach@0 288
msach@0 289 /* this is run after the VMS is set up*/
msach@0 290 void raytrace(void *pixels, VirtProcr *VProc)
msach@0 291 {
msach@0 292 int i;
msach@0 293 double sl, sl_per_procr;
msach@0 294
msach@0 295 /* initialize the random number tables for the jitter */
msach@0 296 for(i=0; i<NRAN; i++) urand[i].x = (double)rand() / RAND_MAX - 0.5;
msach@0 297 for(i=0; i<NRAN; i++) urand[i].y = (double)rand() / RAND_MAX - 0.5;
msach@0 298 for(i=0; i<NRAN; i++) irand[i] = (int)(NRAN * ((double)rand() / RAND_MAX));
msach@0 299
msach@0 300 if(thread_num > yres) {
msach@0 301 fprintf(stderr, "more threads than scanlines specified, reducing number of threads to %d\n", yres);
msach@0 302 thread_num = yres;
msach@0 303 }
msach@0 304
msach@0 305
msach@0 306 if(!(procrs = SSR__malloc_to(thread_num * sizeof(procr_data), VProc))) {
msach@0 307 perror("failed to allocate thread table");
msach@0 308 exit(EXIT_FAILURE);
msach@0 309 }
msach@0 310
msach@0 311 sl = 0.0;
msach@0 312 sl_per_procr = (double)yres / (double)thread_num;
msach@0 313 for(i=0; i<thread_num; i++) {
msach@0 314 procrs[i].sl_start = (int)sl;
msach@0 315 sl += sl_per_procr;
msach@0 316 procrs[i].sl_count = (int)sl - procrs[i].sl_start;
msach@0 317 procrs[i].pixels = (uint32_t*)pixels;
msach@0 318 procrs[i].parentVP = VProc;
msach@0 319
msach@0 320 procrs[i].VP = SSR__create_procr_with((VirtProcrFnPtr)thread_func,
msach@0 321 (void*)(&procrs[i]), VProc);
msach@0 322 }
msach@0 323
msach@0 324 procrs[thread_num - 1].sl_count = yres - procrs[thread_num - 1].sl_start;
msach@0 325
msach@0 326 fprintf(stderr, VER_STR, VER_MAJOR, VER_MINOR);
msach@0 327
msach@0 328 // start worker threads
msach@0 329 //printf("start of worker thread (%d)\n", VProc->procrID);
msach@0 330 start_time = get_msec();
msach@0 331 for(i=0; i<thread_num; i++)
msach@0 332 SSR__send_of_type_to(VProc, NULL, WORK_START, procrs[i].VP);
msach@0 333
msach@0 334 //printf("wait for worker (%d)\n", VProc->procrID);
msach@0 335 for(i=0; i<thread_num; i++)
msach@0 336 SSR__receive_type_to(WORK_END, VProc);
msach@0 337
msach@0 338 rend_time = get_msec() - start_time;
msach@0 339
msach@0 340 SSR__free(procrs,VProc);
msach@0 341 SSR__dissipate_procr(VProc);
msach@0 342 }
msach@0 343
msach@0 344 /* render a frame of xsz/ysz dimensions into the provided framebuffer */
msach@0 345 void render_scanline(int xsz, int ysz, int sl, uint32_t *fb, int samples) {
msach@0 346 int i, s;
msach@0 347 double rcp_samples = 1.0 / (double)samples;
msach@0 348
msach@0 349 for(i=0; i<xsz; i++) {
msach@0 350 double r, g, b;
msach@0 351 r = g = b = 0.0;
msach@0 352
msach@0 353 for(s=0; s<samples; s++) {
msach@0 354 struct vec3 col = trace(get_primary_ray(i, sl, s), 0);
msach@0 355 r += col.x;
msach@0 356 g += col.y;
msach@0 357 b += col.z;
msach@0 358 }
msach@0 359
msach@0 360 r = r * rcp_samples;
msach@0 361 g = g * rcp_samples;
msach@0 362 b = b * rcp_samples;
msach@0 363
msach@0 364 fb[sl * xsz + i] = ((uint32_t)(MIN(r, 1.0) * 255.0) & 0xff) << RSHIFT |
msach@0 365 ((uint32_t)(MIN(g, 1.0) * 255.0) & 0xff) << GSHIFT |
msach@0 366 ((uint32_t)(MIN(b, 1.0) * 255.0) & 0xff) << BSHIFT;
msach@0 367 }
msach@0 368 }
msach@0 369
msach@0 370 /* trace a ray throught the scene recursively (the recursion happens through
msach@0 371 * shade() to calculate reflection rays if necessary).
msach@0 372 */
msach@0 373 struct vec3 trace(struct ray ray, int depth) {
msach@0 374 struct vec3 col;
msach@0 375 struct spoint sp, nearest_sp;
msach@0 376 struct sphere *nearest_obj = 0;
msach@0 377 struct sphere *iter = obj_list->next;
msach@0 378
msach@0 379 /* if we reached the recursion limit, bail out */
msach@0 380 if(depth >= MAX_RAY_DEPTH) {
msach@0 381 col.x = col.y = col.z = 0.0;
msach@0 382 return col;
msach@0 383 }
msach@0 384
msach@0 385 /* find the nearest intersection ... */
msach@0 386 while(iter) {
msach@0 387 if(ray_sphere(iter, ray, &sp)) {
msach@0 388 if(!nearest_obj || sp.dist < nearest_sp.dist) {
msach@0 389 nearest_obj = iter;
msach@0 390 nearest_sp = sp;
msach@0 391 }
msach@0 392 }
msach@0 393 iter = iter->next;
msach@0 394 }
msach@0 395
msach@0 396 /* and perform shading calculations as needed by calling shade() */
msach@0 397 if(nearest_obj) {
msach@0 398 col = shade(nearest_obj, &nearest_sp, depth);
msach@0 399 } else {
msach@0 400 col.x = col.y = col.z = 0.0;
msach@0 401 }
msach@0 402
msach@0 403 return col;
msach@0 404 }
msach@0 405
msach@0 406 /* Calculates direct illumination with the phong reflectance model.
msach@0 407 * Also handles reflections by calling trace again, if necessary.
msach@0 408 */
msach@0 409 struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth) {
msach@0 410 int i;
msach@0 411 struct vec3 col = {0, 0, 0};
msach@0 412
msach@0 413 /* for all lights ... */
msach@0 414 for(i=0; i<lnum; i++) {
msach@0 415 double ispec, idiff;
msach@0 416 struct vec3 ldir;
msach@0 417 struct ray shadow_ray;
msach@0 418 struct sphere *iter = obj_list->next;
msach@0 419 int in_shadow = 0;
msach@0 420
msach@0 421 ldir.x = lights[i].x - sp->pos.x;
msach@0 422 ldir.y = lights[i].y - sp->pos.y;
msach@0 423 ldir.z = lights[i].z - sp->pos.z;
msach@0 424
msach@0 425 shadow_ray.orig = sp->pos;
msach@0 426 shadow_ray.dir = ldir;
msach@0 427
msach@0 428 /* shoot shadow rays to determine if we have a line of sight with the light */
msach@0 429 while(iter) {
msach@0 430 if(ray_sphere(iter, shadow_ray, 0)) {
msach@0 431 in_shadow = 1;
msach@0 432 break;
msach@0 433 }
msach@0 434 iter = iter->next;
msach@0 435 }
msach@0 436
msach@0 437 /* and if we're not in shadow, calculate direct illumination with the phong model. */
msach@0 438 if(!in_shadow) {
msach@0 439 NORMALIZE(ldir);
msach@0 440
msach@0 441 idiff = MAX(DOT(sp->normal, ldir), 0.0);
msach@0 442 ispec = obj->mat.spow > 0.0 ? pow(MAX(DOT(sp->vref, ldir), 0.0), obj->mat.spow) : 0.0;
msach@0 443
msach@0 444 col.x += idiff * obj->mat.col.x + ispec;
msach@0 445 col.y += idiff * obj->mat.col.y + ispec;
msach@0 446 col.z += idiff * obj->mat.col.z + ispec;
msach@0 447 }
msach@0 448 }
msach@0 449
msach@0 450 /* Also, if the object is reflective, spawn a reflection ray, and call trace()
msach@0 451 * to calculate the light arriving from the mirror direction.
msach@0 452 */
msach@0 453 if(obj->mat.refl > 0.0) {
msach@0 454 struct ray ray;
msach@0 455 struct vec3 rcol;
msach@0 456
msach@0 457 ray.orig = sp->pos;
msach@0 458 ray.dir = sp->vref;
msach@0 459 ray.dir.x *= RAY_MAG;
msach@0 460 ray.dir.y *= RAY_MAG;
msach@0 461 ray.dir.z *= RAY_MAG;
msach@0 462
msach@0 463 rcol = trace(ray, depth + 1);
msach@0 464 col.x += rcol.x * obj->mat.refl;
msach@0 465 col.y += rcol.y * obj->mat.refl;
msach@0 466 col.z += rcol.z * obj->mat.refl;
msach@0 467 }
msach@0 468
msach@0 469 return col;
msach@0 470 }
msach@0 471
msach@0 472 /* calculate reflection vector */
msach@0 473 struct vec3 reflect(struct vec3 v, struct vec3 n) {
msach@0 474 struct vec3 res;
msach@0 475 double dot = v.x * n.x + v.y * n.y + v.z * n.z;
msach@0 476 res.x = -(2.0 * dot * n.x - v.x);
msach@0 477 res.y = -(2.0 * dot * n.y - v.y);
msach@0 478 res.z = -(2.0 * dot * n.z - v.z);
msach@0 479 return res;
msach@0 480 }
msach@0 481
msach@0 482 struct vec3 cross_product(struct vec3 v1, struct vec3 v2) {
msach@0 483 struct vec3 res;
msach@0 484 res.x = v1.y * v2.z - v1.z * v2.y;
msach@0 485 res.y = v1.z * v2.x - v1.x * v2.z;
msach@0 486 res.z = v1.x * v2.y - v1.y * v2.x;
msach@0 487 return res;
msach@0 488 }
msach@0 489
msach@0 490 /* determine the primary ray corresponding to the specified pixel (x, y) */
msach@0 491 struct ray get_primary_ray(int x, int y, int sample) {
msach@0 492 struct ray ray;
msach@0 493 float m[3][3];
msach@0 494 struct vec3 i, j = {0, 1, 0}, k, dir, orig, foo;
msach@0 495
msach@0 496 k.x = cam.targ.x - cam.pos.x;
msach@0 497 k.y = cam.targ.y - cam.pos.y;
msach@0 498 k.z = cam.targ.z - cam.pos.z;
msach@0 499 NORMALIZE(k);
msach@0 500
msach@0 501 i = cross_product(j, k);
msach@0 502 j = cross_product(k, i);
msach@0 503 m[0][0] = i.x; m[0][1] = j.x; m[0][2] = k.x;
msach@0 504 m[1][0] = i.y; m[1][1] = j.y; m[1][2] = k.y;
msach@0 505 m[2][0] = i.z; m[2][1] = j.z; m[2][2] = k.z;
msach@0 506
msach@0 507 ray.orig.x = ray.orig.y = ray.orig.z = 0.0;
msach@0 508 ray.dir = get_sample_pos(x, y, sample);
msach@0 509 ray.dir.z = 1.0 / HALF_FOV;
msach@0 510 ray.dir.x *= RAY_MAG;
msach@0 511 ray.dir.y *= RAY_MAG;
msach@0 512 ray.dir.z *= RAY_MAG;
msach@0 513
msach@0 514 dir.x = ray.dir.x + ray.orig.x;
msach@0 515 dir.y = ray.dir.y + ray.orig.y;
msach@0 516 dir.z = ray.dir.z + ray.orig.z;
msach@0 517 foo.x = dir.x * m[0][0] + dir.y * m[0][1] + dir.z * m[0][2];
msach@0 518 foo.y = dir.x * m[1][0] + dir.y * m[1][1] + dir.z * m[1][2];
msach@0 519 foo.z = dir.x * m[2][0] + dir.y * m[2][1] + dir.z * m[2][2];
msach@0 520
msach@0 521 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 522 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 523 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 524
msach@0 525 ray.orig = orig;
msach@0 526 ray.dir.x = foo.x + orig.x;
msach@0 527 ray.dir.y = foo.y + orig.y;
msach@0 528 ray.dir.z = foo.z + orig.z;
msach@0 529
msach@0 530 return ray;
msach@0 531 }
msach@0 532
msach@0 533
msach@0 534 struct vec3 get_sample_pos(int x, int y, int sample) {
msach@0 535 struct vec3 pt;
msach@0 536 static double sf = 0.0;
msach@0 537
msach@0 538 if(sf == 0.0) {
msach@0 539 sf = 1.5 / (double)xres;
msach@0 540 }
msach@0 541
msach@0 542 pt.x = ((double)x / (double)xres) - 0.5;
msach@0 543 pt.y = -(((double)y / (double)yres) - 0.65) / aspect;
msach@0 544
msach@0 545 if(sample) {
msach@0 546 struct vec3 jt = jitter(x, y, sample);
msach@0 547 pt.x += jt.x * sf;
msach@0 548 pt.y += jt.y * sf / aspect;
msach@0 549 }
msach@0 550 return pt;
msach@0 551 }
msach@0 552
msach@0 553 /* jitter function taken from Graphics Gems I. */
msach@0 554 struct vec3 jitter(int x, int y, int s) {
msach@0 555 struct vec3 pt;
msach@0 556 pt.x = urand[(x + (y << 2) + irand[(x + s) & MASK]) & MASK].x;
msach@0 557 pt.y = urand[(y + (x << 2) + irand[(y + s) & MASK]) & MASK].y;
msach@0 558 return pt;
msach@0 559 }
msach@0 560
msach@0 561 /* Calculate ray-sphere intersection, and return {1, 0} to signify hit or no hit.
msach@0 562 * Also the surface point parameters like position, normal, etc are returned through
msach@0 563 * the sp pointer if it is not NULL.
msach@0 564 */
msach@0 565 int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp) {
msach@0 566 double a, b, c, d, sqrt_d, t1, t2;
msach@0 567
msach@0 568 a = SQ(ray.dir.x) + SQ(ray.dir.y) + SQ(ray.dir.z);
msach@0 569 b = 2.0 * ray.dir.x * (ray.orig.x - sph->pos.x) +
msach@0 570 2.0 * ray.dir.y * (ray.orig.y - sph->pos.y) +
msach@0 571 2.0 * ray.dir.z * (ray.orig.z - sph->pos.z);
msach@0 572 c = SQ(sph->pos.x) + SQ(sph->pos.y) + SQ(sph->pos.z) +
msach@0 573 SQ(ray.orig.x) + SQ(ray.orig.y) + SQ(ray.orig.z) +
msach@0 574 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 575
msach@0 576 if((d = SQ(b) - 4.0 * a * c) < 0.0) return 0;
msach@0 577
msach@0 578 sqrt_d = sqrt(d);
msach@0 579 t1 = (-b + sqrt_d) / (2.0 * a);
msach@0 580 t2 = (-b - sqrt_d) / (2.0 * a);
msach@0 581
msach@0 582 if((t1 < ERR_MARGIN && t2 < ERR_MARGIN) || (t1 > 1.0 && t2 > 1.0)) return 0;
msach@0 583
msach@0 584 if(sp) {
msach@0 585 if(t1 < ERR_MARGIN) t1 = t2;
msach@0 586 if(t2 < ERR_MARGIN) t2 = t1;
msach@0 587 sp->dist = t1 < t2 ? t1 : t2;
msach@0 588
msach@0 589 sp->pos.x = ray.orig.x + ray.dir.x * sp->dist;
msach@0 590 sp->pos.y = ray.orig.y + ray.dir.y * sp->dist;
msach@0 591 sp->pos.z = ray.orig.z + ray.dir.z * sp->dist;
msach@0 592
msach@0 593 sp->normal.x = (sp->pos.x - sph->pos.x) / sph->rad;
msach@0 594 sp->normal.y = (sp->pos.y - sph->pos.y) / sph->rad;
msach@0 595 sp->normal.z = (sp->pos.z - sph->pos.z) / sph->rad;
msach@0 596
msach@0 597 sp->vref = reflect(ray.dir, sp->normal);
msach@0 598 NORMALIZE(sp->vref);
msach@0 599 }
msach@0 600 return 1;
msach@0 601 }
msach@0 602
msach@0 603 /* Load the scene from an extremely simple scene description file */
msach@0 604 #define DELIM " \t\n"
msach@0 605 void load_scene(FILE *fp) {
msach@0 606 char line[256], *ptr, type;
msach@0 607
msach@0 608 obj_list = malloc(sizeof(struct sphere));
msach@0 609 obj_list->next = 0;
msach@0 610
msach@0 611 while((ptr = fgets(line, 256, fp))) {
msach@0 612 int i;
msach@0 613 struct vec3 pos, col;
msach@0 614 double rad, spow, refl;
msach@0 615
msach@0 616 while(*ptr == ' ' || *ptr == '\t') ptr++;
msach@0 617 if(*ptr == '#' || *ptr == '\n') continue;
msach@0 618
msach@0 619 if(!(ptr = strtok(line, DELIM))) continue;
msach@0 620 type = *ptr;
msach@0 621
msach@0 622 for(i=0; i<3; i++) {
msach@0 623 if(!(ptr = strtok(0, DELIM))) break;
msach@0 624 *((double*)&pos.x + i) = atof(ptr);
msach@0 625 }
msach@0 626
msach@0 627 if(type == 'l') {
msach@0 628 lights[lnum++] = pos;
msach@0 629 continue;
msach@0 630 }
msach@0 631
msach@0 632 if(!(ptr = strtok(0, DELIM))) continue;
msach@0 633 rad = atof(ptr);
msach@0 634
msach@0 635 for(i=0; i<3; i++) {
msach@0 636 if(!(ptr = strtok(0, DELIM))) break;
msach@0 637 *((double*)&col.x + i) = atof(ptr);
msach@0 638 }
msach@0 639
msach@0 640 if(type == 'c') {
msach@0 641 cam.pos = pos;
msach@0 642 cam.targ = col;
msach@0 643 cam.fov = rad;
msach@0 644 continue;
msach@0 645 }
msach@0 646
msach@0 647 if(!(ptr = strtok(0, DELIM))) continue;
msach@0 648 spow = atof(ptr);
msach@0 649
msach@0 650 if(!(ptr = strtok(0, DELIM))) continue;
msach@0 651 refl = atof(ptr);
msach@0 652
msach@0 653 if(type == 's') {
msach@0 654 struct sphere *sph = malloc(sizeof *sph);
msach@0 655 sph->next = obj_list->next;
msach@0 656 obj_list->next = sph;
msach@0 657
msach@0 658 sph->pos = pos;
msach@0 659 sph->rad = rad;
msach@0 660 sph->mat.col = col;
msach@0 661 sph->mat.spow = spow;
msach@0 662 sph->mat.refl = refl;
msach@0 663 } else {
msach@0 664 fprintf(stderr, "unknown type: %c\n", type);
msach@0 665 }
msach@0 666 }
msach@0 667 }
msach@0 668
msach@0 669
msach@0 670 /* provide a millisecond-resolution timer for each system */
msach@0 671 #if defined(unix) || defined(__unix__)
msach@0 672 #include <time.h>
msach@0 673 #include <sys/time.h>
msach@0 674 unsigned long get_msec(void) {
msach@0 675 static struct timeval timeval, first_timeval;
msach@0 676
msach@0 677 gettimeofday(&timeval, 0);
msach@0 678 if(first_timeval.tv_sec == 0) {
msach@0 679 first_timeval = timeval;
msach@0 680 return 0;
msach@0 681 }
msach@0 682 return (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000;
msach@0 683 }
msach@0 684 #elif defined(__WIN32__) || defined(WIN32)
msach@0 685 #include <windows.h>
msach@0 686 unsigned long get_msec(void) {
msach@0 687 return GetTickCount();
msach@0 688 }
msach@0 689 #else
msach@0 690 #error "I don't know how to measure time on your platform"
msach@0 691 #endif
msach@0 692
msach@0 693 void thread_func(void *tdata, VirtProcr *VProc) {
msach@0 694 int i;
msach@0 695 procr_data *td = (procr_data*)tdata;
msach@0 696
msach@0 697 SSR__receive_type_to(WORK_START, VProc);
msach@0 698
msach@0 699 for(i=0; i<td->sl_count; i++) {
msach@0 700 render_scanline(xres, yres, i + td->sl_start, td->pixels, rays_per_pixel);
msach@0 701 }
msach@0 702
msach@0 703 SSR__send_of_type_to(VProc, NULL, WORK_END, td->parentVP);
msach@0 704
msach@0 705 SSR__dissipate_procr(VProc);
msach@0 706 }