1 /* c-ray-f - a simple raytracing filter.
2 * Copyright (C) 2006 John Tsiombikas <nuclear@siggraph.org>
3 *
4 * You are free to use, modify and redistribute this program under the
5 * terms of the GNU General Public License v2 or (at your option) later.
6 * see "http://www.gnu.org/licenses/gpl.txt" for details.
7 * ---------------------------------------------------------------------
8 * Usage:
9 * compile: cc -o c-ray-f c-ray-f.c -lm
10 * run: cat scene | ./c-ray-f >foo.ppm
11 * enjoy: display foo.ppm (with imagemagick)
12 * or: imgview foo.ppm (on IRIX)
13 * ---------------------------------------------------------------------
14 * Scene file format:
15 * # sphere (many)
16 * s x y z rad r g b shininess reflectivity
17 * # light (many)
18 * l x y z
19 * # camera (one)
20 * c x y z fov tx ty tz
21 * ---------------------------------------------------------------------
22 */
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <math.h>
27 #include <ctype.h>
28 #include <errno.h>
29
30 /* find the appropriate way to define explicitly sized types */
31 #if (__STDC_VERSION__ >= 199900) || defined(__GLIBC__) /* C99 or GNU libc */
32 #include <stdint.h>
33 #elif defined(__unix__) || defined(unix)
34 #include <sys/types.h>
35 #elif defined(_MSC_VER) /* the nameless one */
36 typedef unsigned __int8 uint8_t;
37 typedef unsigned __int32 uint32_t;
38 #endif
39
40 struct vec3 {
41 double x, y, z;
42 };
43
44 struct ray {
45 struct vec3 orig, dir;
46 };
47
48 struct material {
49 struct vec3 col; /* color */
50 double spow; /* specular power */
51 double refl; /* reflection intensity */
52 };
53
54 struct sphere {
55 struct vec3 pos;
56 double rad;
57 struct material mat;
58 struct sphere *next;
59 };
60
61 struct spoint {
62 struct vec3 pos, normal, vref; /* position, normal and view reflection */
63 double dist; /* parametric distance of intersection along the ray */
64 };
65
66 struct camera {
67 struct vec3 pos, targ;
68 double fov;
69 };
70
71 void render(int xsz, int ysz, uint32_t *fb, int samples);
72 struct vec3 trace(struct ray ray, int depth);
73 struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth);
74 struct vec3 reflect(struct vec3 v, struct vec3 n);
75 struct vec3 cross_product(struct vec3 v1, struct vec3 v2);
76 struct ray get_primary_ray(int x, int y, int sample);
77 struct vec3 get_sample_pos(int x, int y, int sample);
78 struct vec3 jitter(int x, int y, int s);
79 int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp);
80 void load_scene(FILE *fp);
81 unsigned long get_msec(void);
82
83 #define MAX_LIGHTS 16 /* maximum number of lights */
84 #define RAY_MAG 1000.0 /* trace rays of this magnitude */
85 #define MAX_RAY_DEPTH 5 /* raytrace recursion limit */
86 #define FOV 0.78539816 /* field of view in rads (pi/4) */
87 #define HALF_FOV (FOV * 0.5)
88 #define ERR_MARGIN 1e-6 /* an arbitrary error margin to avoid surface acne */
89
90 /* bit-shift amount for packing each color into a 32bit uint */
91 #ifdef LITTLE_ENDIAN
92 #define RSHIFT 16
93 #define BSHIFT 0
94 #else /* big endian */
95 #define RSHIFT 0
96 #define BSHIFT 16
97 #endif /* endianness */
98 #define GSHIFT 8 /* this is the same in both byte orders */
99
100 /* some helpful macros... */
101 #define SQ(x) ((x) * (x))
102 #define MAX(a, b) ((a) > (b) ? (a) : (b))
103 #define MIN(a, b) ((a) < (b) ? (a) : (b))
104 #define DOT(a, b) ((a).x * (b).x + (a).y * (b).y + (a).z * (b).z)
105 #define NORMALIZE(a) do {\
106 double len = sqrt(DOT(a, a));\
107 (a).x /= len; (a).y /= len; (a).z /= len;\
108 } while(0);
109
110 /* global state */
111 int xres = 800;
112 int yres = 600;
113 double aspect = 1.333333;
114 struct sphere *obj_list;
115 struct vec3 lights[MAX_LIGHTS];
116 int lnum = 0;
117 struct camera cam;
118
119 #define NRAN 1024
120 #define MASK (NRAN - 1)
121 struct vec3 urand[NRAN];
122 int irand[NRAN];
123
124 const char *usage = {
125 "Usage: c-ray-f [options]\n"
126 " Reads a scene file from stdin, writes the image to stdout, and stats to stderr.\n\n"
127 "Options:\n"
128 " -s WxH where W is the width and H the height of the image\n"
129 " -r <rays> shoot <rays> rays per pixel (antialiasing)\n"
130 " -i <file> read from <file> instead of stdin\n"
131 " -o <file> write to <file> instead of stdout\n"
132 " -h this help screen\n\n"
133 };
134
135
136
main(int argc,char ** argv)137 int main(int argc, char **argv) {
138 int i;
139 unsigned long rend_time, start_time;
140 uint32_t *pixels;
141 int rays_per_pixel = 1;
142 FILE *infile = stdin, *outfile = stdout;
143
144 for(i=1; i<argc; i++) {
145 if(argv[i][0] == '-' && argv[i][2] == 0) {
146 char *sep;
147 switch(argv[i][1]) {
148 case 's':
149 if(!isdigit(argv[++i][0]) || !(sep = strchr(argv[i], 'x')) || !isdigit(*(sep + 1))) {
150 fputs("-s must be followed by something like \"640x480\"\n", stderr);
151 return EXIT_FAILURE;
152 }
153 xres = atoi(argv[i]);
154 yres = atoi(sep + 1);
155 aspect = (double)xres / (double)yres;
156 break;
157
158 case 'i':
159 if(!(infile = fopen(argv[++i], "r"))) {
160 fprintf(stderr, "failed to open input file %s: %s\n", argv[i], strerror(errno));
161 return EXIT_FAILURE;
162 }
163 break;
164
165 case 'o':
166 if(!(outfile = fopen(argv[++i], "w"))) {
167 fprintf(stderr, "failed to open output file %s: %s\n", argv[i], strerror(errno));
168 return EXIT_FAILURE;
169 }
170 break;
171
172 case 'r':
173 if(!isdigit(argv[++i][0])) {
174 fputs("-r must be followed by a number (rays per pixel)\n", stderr);
175 return EXIT_FAILURE;
176 }
177 rays_per_pixel = atoi(argv[i]);
178 break;
179
180 case 'h':
181 fputs(usage, stdout);
182 return 0;
183
184 default:
185 fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
186 fputs(usage, stderr);
187 return EXIT_FAILURE;
188 }
189 } else {
190 fprintf(stderr, "unrecognized argument: %s\n", argv[i]);
191 fputs(usage, stderr);
192 return EXIT_FAILURE;
193 }
194 }
195
196 if(!(pixels = malloc(xres * yres * sizeof *pixels))) {
197 perror("pixel buffer allocation failed");
198 return EXIT_FAILURE;
199 }
200 load_scene(infile);
201
202 /* initialize the random number tables for the jitter */
203 for(i=0; i<NRAN; i++) urand[i].x = (double)rand() / RAND_MAX - 0.5;
204 for(i=0; i<NRAN; i++) urand[i].y = (double)rand() / RAND_MAX - 0.5;
205 for(i=0; i<NRAN; i++) irand[i] = (int)(NRAN * ((double)rand() / RAND_MAX));
206
207 start_time = get_msec();
208 render(xres, yres, pixels, rays_per_pixel);
209 rend_time = get_msec() - start_time;
210
211 /* output statistics to stderr */
212 fprintf(stderr, "Rendering took: %lu seconds (%lu milliseconds)\n", rend_time / 1000, rend_time);
213
214 /* output the image */
215 fprintf(outfile, "P6\n%d %d\n255\n", xres, yres);
216 for(i=0; i<xres * yres; i++) {
217 fputc((pixels[i] >> RSHIFT) & 0xff, outfile);
218 fputc((pixels[i] >> GSHIFT) & 0xff, outfile);
219 fputc((pixels[i] >> BSHIFT) & 0xff, outfile);
220 }
221 fflush(outfile);
222
223 if(infile != stdin) fclose(infile);
224 if(outfile != stdout) fclose(outfile);
225 return 0;
226 }
227
228 /* render a frame of xsz/ysz dimensions into the provided framebuffer */
render(int xsz,int ysz,uint32_t * fb,int samples)229 void render(int xsz, int ysz, uint32_t *fb, int samples) {
230 int i, j, s;
231 double rcp_samples = 1.0 / (double)samples;
232
233 /* for each subpixel, trace a ray through the scene, accumulate the
234 * colors of the subpixels of each pixel, then pack the color and
235 * put it into the framebuffer.
236 * XXX: assumes contiguous scanlines with NO padding, and 32bit pixels.
237 */
238 for(j=0; j<ysz; j++) {
239 for(i=0; i<xsz; i++) {
240 double r, g, b;
241 r = g = b = 0.0;
242
243 for(s=0; s<samples; s++) {
244 struct vec3 col = trace(get_primary_ray(i, j, s), 0);
245 r += col.x;
246 g += col.y;
247 b += col.z;
248 }
249
250 r = r * rcp_samples;
251 g = g * rcp_samples;
252 b = b * rcp_samples;
253
254 *fb++ = ((uint32_t)(MIN(r, 1.0) * 255.0) & 0xff) << RSHIFT |
255 ((uint32_t)(MIN(g, 1.0) * 255.0) & 0xff) << GSHIFT |
256 ((uint32_t)(MIN(b, 1.0) * 255.0) & 0xff) << BSHIFT;
257 }
258 }
259 }
260
261 /* trace a ray through the scene recursively (the recursion happens through
262 * shade() to calculate reflection rays if necessary).
263 */
trace(struct ray ray,int depth)264 struct vec3 trace(struct ray ray, int depth) {
265 struct vec3 col;
266 struct spoint sp, nearest_sp;
267 struct sphere *nearest_obj = 0;
268 struct sphere *iter = obj_list->next;
269
270 /* if we reached the recursion limit, bail out */
271 if(depth >= MAX_RAY_DEPTH) {
272 col.x = col.y = col.z = 0.0;
273 return col;
274 }
275
276 /* find the nearest intersection ... */
277 while(iter) {
278 if(ray_sphere(iter, ray, &sp)) {
279 if(!nearest_obj || sp.dist < nearest_sp.dist) {
280 nearest_obj = iter;
281 nearest_sp = sp;
282 }
283 }
284 iter = iter->next;
285 }
286
287 /* and perform shading calculations as needed by calling shade() */
288 if(nearest_obj) {
289 col = shade(nearest_obj, &nearest_sp, depth);
290 } else {
291 col.x = col.y = col.z = 0.0;
292 }
293
294 return col;
295 }
296
297 /* Calculates direct illumination with the phong reflectance model.
298 * Also handles reflections by calling trace again, if necessary.
299 */
shade(struct sphere * obj,struct spoint * sp,int depth)300 struct vec3 shade(struct sphere *obj, struct spoint *sp, int depth) {
301 int i;
302 struct vec3 col = {0, 0, 0};
303
304 /* for all lights ... */
305 for(i=0; i<lnum; i++) {
306 double ispec, idiff;
307 struct vec3 ldir;
308 struct ray shadow_ray;
309 struct sphere *iter = obj_list->next;
310 int in_shadow = 0;
311
312 ldir.x = lights[i].x - sp->pos.x;
313 ldir.y = lights[i].y - sp->pos.y;
314 ldir.z = lights[i].z - sp->pos.z;
315
316 shadow_ray.orig = sp->pos;
317 shadow_ray.dir = ldir;
318
319 /* shoot shadow rays to determine if we have a line of sight with the light */
320 while(iter) {
321 if(ray_sphere(iter, shadow_ray, 0)) {
322 in_shadow = 1;
323 break;
324 }
325 iter = iter->next;
326 }
327
328 /* and if we're not in shadow, calculate direct illumination with the phong model. */
329 if(!in_shadow) {
330 NORMALIZE(ldir);
331
332 idiff = MAX(DOT(sp->normal, ldir), 0.0);
333 ispec = obj->mat.spow > 0.0 ? pow(MAX(DOT(sp->vref, ldir), 0.0), obj->mat.spow) : 0.0;
334
335 col.x += idiff * obj->mat.col.x + ispec;
336 col.y += idiff * obj->mat.col.y + ispec;
337 col.z += idiff * obj->mat.col.z + ispec;
338 }
339 }
340
341 /* Also, if the object is reflective, spawn a reflection ray, and call trace()
342 * to calculate the light arriving from the mirror direction.
343 */
344 if(obj->mat.refl > 0.0) {
345 struct ray ray;
346 struct vec3 rcol;
347
348 ray.orig = sp->pos;
349 ray.dir = sp->vref;
350 ray.dir.x *= RAY_MAG;
351 ray.dir.y *= RAY_MAG;
352 ray.dir.z *= RAY_MAG;
353
354 rcol = trace(ray, depth + 1);
355 col.x += rcol.x * obj->mat.refl;
356 col.y += rcol.y * obj->mat.refl;
357 col.z += rcol.z * obj->mat.refl;
358 }
359
360 return col;
361 }
362
363 /* calculate reflection vector */
reflect(struct vec3 v,struct vec3 n)364 struct vec3 reflect(struct vec3 v, struct vec3 n) {
365 struct vec3 res;
366 double dot = v.x * n.x + v.y * n.y + v.z * n.z;
367 res.x = -(2.0 * dot * n.x - v.x);
368 res.y = -(2.0 * dot * n.y - v.y);
369 res.z = -(2.0 * dot * n.z - v.z);
370 return res;
371 }
372
cross_product(struct vec3 v1,struct vec3 v2)373 struct vec3 cross_product(struct vec3 v1, struct vec3 v2) {
374 struct vec3 res;
375 res.x = v1.y * v2.z - v1.z * v2.y;
376 res.y = v1.z * v2.x - v1.x * v2.z;
377 res.z = v1.x * v2.y - v1.y * v2.x;
378 return res;
379 }
380
381 /* determine the primary ray corresponding to the specified pixel (x, y) */
get_primary_ray(int x,int y,int sample)382 struct ray get_primary_ray(int x, int y, int sample) {
383 struct ray ray;
384 float m[3][3];
385 struct vec3 i, j = {0, 1, 0}, k, dir, orig, foo;
386
387 k.x = cam.targ.x - cam.pos.x;
388 k.y = cam.targ.y - cam.pos.y;
389 k.z = cam.targ.z - cam.pos.z;
390 NORMALIZE(k);
391
392 i = cross_product(j, k);
393 j = cross_product(k, i);
394 m[0][0] = i.x; m[0][1] = j.x; m[0][2] = k.x;
395 m[1][0] = i.y; m[1][1] = j.y; m[1][2] = k.y;
396 m[2][0] = i.z; m[2][1] = j.z; m[2][2] = k.z;
397
398 ray.orig.x = ray.orig.y = ray.orig.z = 0.0;
399 ray.dir = get_sample_pos(x, y, sample);
400 ray.dir.z = 1.0 / HALF_FOV;
401 ray.dir.x *= RAY_MAG;
402 ray.dir.y *= RAY_MAG;
403 ray.dir.z *= RAY_MAG;
404
405 dir.x = ray.dir.x + ray.orig.x;
406 dir.y = ray.dir.y + ray.orig.y;
407 dir.z = ray.dir.z + ray.orig.z;
408 foo.x = dir.x * m[0][0] + dir.y * m[0][1] + dir.z * m[0][2];
409 foo.y = dir.x * m[1][0] + dir.y * m[1][1] + dir.z * m[1][2];
410 foo.z = dir.x * m[2][0] + dir.y * m[2][1] + dir.z * m[2][2];
411
412 orig.x = ray.orig.x * m[0][0] + ray.orig.y * m[0][1] + ray.orig.z * m[0][2] + cam.pos.x;
413 orig.y = ray.orig.x * m[1][0] + ray.orig.y * m[1][1] + ray.orig.z * m[1][2] + cam.pos.y;
414 orig.z = ray.orig.x * m[2][0] + ray.orig.y * m[2][1] + ray.orig.z * m[2][2] + cam.pos.z;
415
416 ray.orig = orig;
417 ray.dir.x = foo.x + orig.x;
418 ray.dir.y = foo.y + orig.y;
419 ray.dir.z = foo.z + orig.z;
420
421 return ray;
422 }
423
424
get_sample_pos(int x,int y,int sample)425 struct vec3 get_sample_pos(int x, int y, int sample) {
426 struct vec3 pt;
427 double xsz = 2.0, ysz = xres / aspect;
428 static double sf = 0.0;
429
430 if(sf == 0.0) {
431 sf = 2.0 / (double)xres;
432 }
433
434 pt.x = ((double)x / (double)xres) - 0.5;
435 pt.y = -(((double)y / (double)yres) - 0.65) / aspect;
436
437 if(sample) {
438 struct vec3 jt = jitter(x, y, sample);
439 pt.x += jt.x * sf;
440 pt.y += jt.y * sf / aspect;
441 }
442 return pt;
443 }
444
445 /* jitter function taken from Graphics Gems I. */
jitter(int x,int y,int s)446 struct vec3 jitter(int x, int y, int s) {
447 struct vec3 pt;
448 pt.x = urand[(x + (y << 2) + irand[(x + s) & MASK]) & MASK].x;
449 pt.y = urand[(y + (x << 2) + irand[(y + s) & MASK]) & MASK].y;
450 return pt;
451 }
452
453 /* Calculate ray-sphere intersection, and return {1, 0} to signify hit or no hit.
454 * Also the surface point parameters like position, normal, etc are returned through
455 * the sp pointer if it is not NULL.
456 */
ray_sphere(const struct sphere * sph,struct ray ray,struct spoint * sp)457 int ray_sphere(const struct sphere *sph, struct ray ray, struct spoint *sp) {
458 double a, b, c, d, sqrt_d, t1, t2;
459
460 a = SQ(ray.dir.x) + SQ(ray.dir.y) + SQ(ray.dir.z);
461 b = 2.0 * ray.dir.x * (ray.orig.x - sph->pos.x) +
462 2.0 * ray.dir.y * (ray.orig.y - sph->pos.y) +
463 2.0 * ray.dir.z * (ray.orig.z - sph->pos.z);
464 c = SQ(sph->pos.x) + SQ(sph->pos.y) + SQ(sph->pos.z) +
465 SQ(ray.orig.x) + SQ(ray.orig.y) + SQ(ray.orig.z) +
466 2.0 * (-sph->pos.x * ray.orig.x - sph->pos.y * ray.orig.y - sph->pos.z * ray.orig.z) - SQ(sph->rad);
467
468 if((d = SQ(b) - 4.0 * a * c) < 0.0) return 0;
469
470 sqrt_d = sqrt(d);
471 t1 = (-b + sqrt_d) / (2.0 * a);
472 t2 = (-b - sqrt_d) / (2.0 * a);
473
474 if((t1 < ERR_MARGIN && t2 < ERR_MARGIN) || (t1 > 1.0 && t2 > 1.0)) return 0;
475
476 if(sp) {
477 if(t1 < ERR_MARGIN) t1 = t2;
478 if(t2 < ERR_MARGIN) t2 = t1;
479 sp->dist = t1 < t2 ? t1 : t2;
480
481 sp->pos.x = ray.orig.x + ray.dir.x * sp->dist;
482 sp->pos.y = ray.orig.y + ray.dir.y * sp->dist;
483 sp->pos.z = ray.orig.z + ray.dir.z * sp->dist;
484
485 sp->normal.x = (sp->pos.x - sph->pos.x) / sph->rad;
486 sp->normal.y = (sp->pos.y - sph->pos.y) / sph->rad;
487 sp->normal.z = (sp->pos.z - sph->pos.z) / sph->rad;
488
489 sp->vref = reflect(ray.dir, sp->normal);
490 NORMALIZE(sp->vref);
491 }
492 return 1;
493 }
494
495 /* Load the scene from an extremely simple scene description file */
496 #define DELIM " \t\n"
load_scene(FILE * fp)497 void load_scene(FILE *fp) {
498 char line[256], *ptr, type;
499
500 obj_list = malloc(sizeof(struct sphere));
501 obj_list->next = 0;
502
503 while((ptr = fgets(line, 256, fp))) {
504 int i;
505 struct vec3 pos, col;
506 double rad, spow, refl;
507
508 while(*ptr == ' ' || *ptr == '\t') ptr++;
509 if(*ptr == '#' || *ptr == '\n') continue;
510
511 if(!(ptr = strtok(line, DELIM))) continue;
512 type = *ptr;
513
514 for(i=0; i<3; i++) {
515 if(!(ptr = strtok(0, DELIM))) break;
516 *((double*)&pos.x + i) = atof(ptr);
517 }
518
519 if(type == 'l') {
520 lights[lnum++] = pos;
521 continue;
522 }
523
524 if(!(ptr = strtok(0, DELIM))) continue;
525 rad = atof(ptr);
526
527 for(i=0; i<3; i++) {
528 if(!(ptr = strtok(0, DELIM))) break;
529 *((double*)&col.x + i) = atof(ptr);
530 }
531
532 if(type == 'c') {
533 cam.pos = pos;
534 cam.targ = col;
535 cam.fov = rad;
536 continue;
537 }
538
539 if(!(ptr = strtok(0, DELIM))) continue;
540 spow = atof(ptr);
541
542 if(!(ptr = strtok(0, DELIM))) continue;
543 refl = atof(ptr);
544
545 if(type == 's') {
546 struct sphere *sph = malloc(sizeof *sph);
547 sph->next = obj_list->next;
548 obj_list->next = sph;
549
550 sph->pos = pos;
551 sph->rad = rad;
552 sph->mat.col = col;
553 sph->mat.spow = spow;
554 sph->mat.refl = refl;
555 } else {
556 fprintf(stderr, "unknown type: %c\n", type);
557 }
558 }
559 }
560
561
562 /* provide a millisecond-resolution timer for each system */
563 #if defined(__unix__) || defined(unix)
564 #include <time.h>
565 #include <sys/time.h>
get_msec(void)566 unsigned long get_msec(void) {
567 static struct timeval timeval, first_timeval;
568
569 gettimeofday(&timeval, 0);
570 if(first_timeval.tv_sec == 0) {
571 first_timeval = timeval;
572 return 0;
573 }
574 return (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000;
575 }
576 #elif defined(__WIN32__) || defined(WIN32)
577 #include <windows.h>
get_msec(void)578 unsigned long get_msec(void) {
579 return GetTickCount();
580 }
581 #else
582 #error "I don't know how to measure time on your platform"
583 #endif
584