examples/readme_loop/ray: alternate version for the READMEHEADmaster

1 files changed, 223 insertions, 0 deletions

diff --git a/examples/readme_loop/ray.c b/examples/readme_loop/ray.c
new file mode 100644
index 0000000..17f81d4
--- /dev/null
+++ b/examples/readme_loop/ray.c
@@ -0,0 +1,223 @@
+#include "ray.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <pthread.h>
+#include <unistd.h>
+
+#include "3dmath.h"
+
+#define LENGTH(array) (sizeof(array) / sizeof(array[0]))
+#define MAX(x, y) (x > y ? x : y)
+#define MIN(x, y) (x < y ? x : y)
+
+#define TAU 6.28318531
+
+#if __GNUC__ >= 3
+#  define unlikely(cond) __builtin_expect ((cond), 0)
+#  define likely(cond)   __builtin_expect ((cond), 1)
+#else
+#  define unlikely(cond) (cond)
+#  define likely(cond) (cond)
+#endif
+
+typedef struct {
+    float position[3];
+    float radius;
+    float diffuse[3];
+    float specular[3];
+    int subtract;
+} Object;
+
+typedef struct {
+    float position[3];
+    float diffuse[3];
+} Light;
+
+typedef struct {
+    pthread_mutex_t mutex;
+
+    int width, height;
+    unsigned char* buffer;
+    long next_line;
+} ThreadArg;
+
+static float* trace_vectors;
+static int trace_vectors_width, trace_vectors_height;
+
+static Object objects[] = {
+    {.position={-1.414, -1, -3}, .radius=1, .diffuse={.8, .0, .8}, .specular={.7, .6, .7}, .subtract=0},
+    {.position={0, 1.414, -3}, .radius=1, .diffuse={.0, .8, .8}, .specular={.6, .7, .7}, .subtract=0},
+    {.position={0, 0, -3}, .radius=1.5, .diffuse={.8, .8, .8}, .specular={.7, .7, .7}, .subtract=1},
+    {.position={1.414, -1, -3}, .radius=1, .diffuse={.8, .8, .0}, .specular={.7, .7, .6}, .subtract=0},
+    {.position={0, 0, -3}, .radius=1.1, .diffuse={.9, .9, .9}, .specular={.9, .9, .9}, .subtract=2}
+};
+static const Light lights[] = {
+    {.position={-3, 3, -4}, .diffuse={0, .6, .6}},
+    {.position={0, 30, -4}, .diffuse={1, 1, 1}}
+};
+static const float ambient[3] = {0.2, 0.1, 0.1};
+
+static void
+trace(const float s[3], const float d[3], float pixel[3], int n) {
+    // Reflections in concave objects can go really deep, so we need to limit
+    // the recursion depth.
+    if (n > 6) return;
+
+    float nearest = HUGE_VAL;
+    int nearest_object = -1;
+    float nearest_y[3];
+    float nearest_r[3];
+
+    for(size_t j = 0; j < LENGTH(objects); ++j) {
+        float r[3], t, y[3];
+
+        if (objects[j].subtract == 1) continue;
+
+        t = sphere_intersect(y, r, s, d, objects[j].position, objects[j].radius, 0);
+
+        if(likely(t <= 0) || t > nearest)
+          continue;
+
+        if (objects[j].subtract == 0) {
+          size_t k;
+          for (k = 0; k < LENGTH(objects); ++k) {
+              if (!objects[k].subtract) continue;
+              if (POW2(y[0] - objects[k].position[0]) + POW2(y[1] - objects[k].position[1]) + POW2(y[2] - objects[k].position[2]) > POW2(objects[k].radius)) continue;
+
+              t = sphere_intersect(y, r, s, d, objects[k].position, objects[k].radius, 1);
+
+              break;
+          }
+
+          if(likely(t <= 0) || t > nearest)
+            continue;
+        }
+
+        nearest = t;
+        nearest_object = j;
+        memcpy(nearest_y, y, sizeof(nearest_y));
+        memcpy(nearest_r, r, sizeof(nearest_y));
+    }
+
+    if (nearest_object == -1) return;
+
+    trace(nearest_y, nearest_r, pixel, n + 1);
+
+    for (int k = 0; k < 3; ++k)
+        pixel[k] = pixel[k] * objects[nearest_object].specular[k] + ambient[k] * objects[nearest_object].diffuse[k];
+
+    for(int m = 0; m < LENGTH(lights); ++m) {
+        float l[3];
+        for(int i = 0; i < 3; ++i)
+            l[i] = lights[m].position[i] - nearest_y[i];
+
+        float lr_dot = dot(l, nearest_r);
+        if (lr_dot <= 0) continue;
+
+        float scale = lr_dot / sqrtf(dot(l, l)) / (1 << n);
+        // The cutoff at 0.05 is for artistic reasons; 0.0 would be more
+        // realistic.
+        if (scale <= 0.05) continue;
+
+        for(int k = 0; k < 3; ++k)
+            pixel[k] += lights[m].diffuse[k] * objects[nearest_object].diffuse[k] * scale;
+    }
+}
+
+static void
+trace_line(int l, int width, unsigned char *buf) {
+    static const float s[3] = {0, 0, 8};
+
+    for(int i = 0; i < width; ++i, buf += 4) {
+        float pixel[3] = { 0, 0, 0 };
+
+        trace(s, &trace_vectors[(l * width + i) * 3], pixel, 1);
+
+        buf[0] = MIN(pixel[0], 1.0f) * 255;
+        buf[1] = MIN(pixel[1], 1.0f) * 255;
+        buf[2] = MIN(pixel[2], 1.0f) * 255;
+    }
+}
+
+static void *
+thread(void *arg) {
+    ThreadArg* thread_arg = arg;
+
+    for (;;) {
+        pthread_mutex_lock(&thread_arg->mutex);
+        if (thread_arg->next_line == thread_arg->height) break;
+        long line = thread_arg->next_line++;
+        pthread_mutex_unlock(&thread_arg->mutex);
+
+        trace_line(line, thread_arg->width, thread_arg->buffer + line * 4 * thread_arg->width);
+    }
+
+    pthread_mutex_unlock(&thread_arg->mutex);
+
+    return NULL;
+}
+
+static void
+initialize_trace_vectors(int width, int height) {
+    trace_vectors = calloc(width * height, 3 * sizeof(float));
+    trace_vectors_width = width;
+    trace_vectors_height = height;
+    for(int y = 0; y < height; ++y) {
+        for(int x = 0; x < width; ++x) {
+          float* d = &trace_vectors[(y * width + x) * 3];
+          d[0] = ((float)x / width - 0.5f) * 0.5f * ((float)width / height);
+          d[1] = ((float)y / height - 0.5f) * 0.5f;
+          d[2] = -1;
+          normalize(d);
+        }
+    }
+}
+
+void
+trace_scene(float time, int width, int height, unsigned char *buf, int threaded) {
+    if (trace_vectors && (trace_vectors_width != width || trace_vectors_height != height)) {
+      free(trace_vectors);
+      trace_vectors = 0;
+    }
+    if (!trace_vectors)
+      initialize_trace_vectors(width, height);
+
+    objects[0].position[0] = (1.5 + 0.35 * sin(1.0 * time + 0.0)) * cos(0.5 * time);
+    objects[0].position[1] = (1.5 + 0.35 * sin(1.0 * time + 2.5)) * sin(0.5 * time);
+    objects[1].position[0] = (1.5 + 0.35 * sin(1.0 * time + 2.0)) * cos(0.5 * time + 1/3. * TAU);
+    objects[1].position[1] = (1.5 + 0.35 * sin(1.0 * time + 1.5)) * sin(0.5 * time + 1/3. * TAU);
+    objects[3].position[0] = (1.5 + 0.35 * sin(1.0 * time + 1.0)) * cos(0.5 * time + 2/3. * TAU);
+    objects[3].position[1] = (1.5 + 0.35 * sin(1.0 * time + 0.5)) * sin(0.5 * time + 2/3. * TAU);
+    objects[2].position[2] = -3 + 0.2 * sin(time * 1.0);
+    memcpy(objects[4].position, objects[2].position, sizeof(objects[4].position));
+
+    if(threaded) {
+        ThreadArg arg;
+        memset(&arg, 0, sizeof(arg));
+        arg.width = width;
+        arg.height = height;
+        pthread_mutex_init(&arg.mutex, NULL);
+        arg.buffer = buf;
+
+        int num_threads = sysconf(_SC_NPROCESSORS_CONF) - 1;
+        pthread_t* threads = NULL;
+        if (num_threads > 0) {
+          threads = calloc(sizeof(*threads), num_threads);
+
+          for (int i = 0; i < num_threads; ++i)
+            pthread_create(&threads[i], NULL, thread, &arg);
+        }
+
+        thread(&arg);
+
+        for(int i = 0; i < num_threads; ++i)
+            pthread_join(threads[i], NULL);
+        free(threads);
+    } else {
+        for(int i = 0; i < height; ++i)
+            trace_line(i, width, buf + i * 4 * width);
+    }
+}