box2.c

#include <assert.h> // so I wont have to do error handling
#include <math.h>
#include <png.h>   // to output frames as png files
#include <stdio.h> // I/O
#include <stdlib.h>
#include <string.h>

#define FRAMES 50

/* this may be a vertex or a vector */
struct vec3D {
  double x, y, z;
};

/* A triangle is composed of 3D vertices */
struct triangle3D {
  struct vec3D vertex[3];
};

/* An object is composed of triangles */
struct object3D {
  int faces;
  struct triangle3D *face;
};

/* A scene is composed of objects */
struct scene3D {
  int objects;
  struct object3D *object;
};

int width, height;
static double dot_product(const struct vec3D a, const struct vec3D b) {
  return a.x * b.x + a.y * b.y + a.z * b.z;
}
static double magnitude(const struct vec3D a) {
  return sqrt(dot_product(a, a));
}
static struct object3D translate(const struct object3D ctx, const double dx,
                                 const double dy, const double dz) {
  struct triangle3D *face = malloc(sizeof(struct triangle3D) * ctx.faces);
  assert(face != NULL);
  for (int i = 0; i < ctx.faces; i++) {
    face[i] =
        (struct triangle3D){.vertex = {[0].x = ctx.face[i].vertex[0].x + dx,
                                       [0].y = ctx.face[i].vertex[0].y + dy,
                                       [0].z = ctx.face[i].vertex[0].z + dz,
                                       [1].x = ctx.face[i].vertex[1].x + dx,
                                       [1].y = ctx.face[i].vertex[1].y + dy,
                                       [1].z = ctx.face[i].vertex[1].z + dz,
                                       [2].x = ctx.face[i].vertex[2].x + dx,
                                       [2].y = ctx.face[i].vertex[2].y + dy,
                                       [2].z = ctx.face[i].vertex[2].z + dz}};
  }
  return (struct object3D){.faces = ctx.faces, .face = face};
}
static struct object3D rotate(const struct object3D ctx, const double dx,
                              const double dy, const double dz) {
  struct triangle3D *face = malloc(sizeof(struct triangle3D) * ctx.faces);
  assert(face != NULL);
  for (int i = 0; i < ctx.faces; i++) {
    face[i] = (struct triangle3D){
        .vertex = {[0].x = cos(dz) * cos(dy) * ctx.face[i].vertex[0].x +
                           (cos(dz) * sin(dy) * sin(dx) - sin(dz) * cos(dx)) *
                               ctx.face[i].vertex[0].y +
                           (cos(dz) * sin(dy) * cos(dx) + sin(dz) * sin(dx)) *
                               ctx.face[i].vertex[0].z,
                   [0].y = sin(dz) * cos(dy) * ctx.face[i].vertex[0].x +
                           (sin(dz) * sin(dy) * sin(dx) + cos(dz) * cos(dx)) *
                               ctx.face[i].vertex[0].y +
                           (sin(dz) * sin(dy) * cos(dx) - cos(dz) * sin(dx)) *
                               ctx.face[i].vertex[0].z,
                   [0].z = -sin(dy) * ctx.face[i].vertex[0].x +
                           cos(dy) * sin(dx) * ctx.face[i].vertex[0].y +
                           cos(dy) * cos(dx) * ctx.face[i].vertex[0].z,

                   [1].x = cos(dz) * cos(dy) * ctx.face[i].vertex[1].x +
                           (cos(dz) * sin(dy) * sin(dx) - sin(dz) * cos(dx)) *
                               ctx.face[i].vertex[1].y +
                           (cos(dz) * sin(dy) * cos(dx) + sin(dz) * sin(dx)) *
                               ctx.face[i].vertex[1].z,
                   [1].y = sin(dz) * cos(dy) * ctx.face[i].vertex[1].x +
                           (sin(dz) * sin(dy) * sin(dx) + cos(dz) * cos(dx)) *
                               ctx.face[i].vertex[1].y +
                           (sin(dz) * sin(dy) * cos(dx) - cos(dz) * sin(dx)) *
                               ctx.face[i].vertex[1].z,
                   [1].z = -sin(dy) * ctx.face[i].vertex[1].x +
                           cos(dy) * sin(dx) * ctx.face[i].vertex[1].y +
                           cos(dy) * cos(dx) * ctx.face[i].vertex[1].z,

                   [2].x = cos(dz) * cos(dy) * ctx.face[i].vertex[2].x +
                           (cos(dz) * sin(dy) * sin(dx) - sin(dz) * cos(dx)) *
                               ctx.face[i].vertex[2].y +
                           (cos(dz) * sin(dy) * cos(dx) + sin(dz) * sin(dx)) *
                               ctx.face[i].vertex[2].z,
                   [2].y = sin(dz) * cos(dy) * ctx.face[i].vertex[2].x +
                           (sin(dz) * sin(dy) * sin(dx) + cos(dz) * cos(dx)) *
                               ctx.face[i].vertex[2].y +
                           (sin(dz) * sin(dy) * cos(dx) - cos(dz) * sin(dx)) *
                               ctx.face[i].vertex[2].z,
                   [2].z = -sin(dy) * ctx.face[i].vertex[2].x +
                           cos(dy) * sin(dx) * ctx.face[i].vertex[2].y +
                           cos(dy) * cos(dx) * ctx.face[i].vertex[2].z}};
  }
  return (struct object3D){.faces = ctx.faces, .face = face};
}
static struct vec3D cross_product(const struct vec3D a, const struct vec3D b) {
  /* designated initializers */
  return (struct vec3D){.x = a.y * b.z - a.z * b.y,
                        .y = a.z * b.x - a.x * b.z,
                        .z = a.x * b.y - a.y * b.x};
}
static struct vec3D normalize(const struct vec3D ctx) {
  double mg = magnitude(ctx);
  return (struct vec3D){.x = ctx.x / mg, .y = ctx.y / mg, .z = ctx.z / mg};
}
int main(int argc, char *argv[]) {
  double Pz = 150.0;
  (void)(width = 1000), height = 720;
  struct triangle3D face[] = {
#include "faces.txt"
  };
  struct object3D box = {.faces = 8, .face = face};
  struct vec3D light_source = {
      .x = 0.0, .y = 0.0, .z = -1.0}; /* magnitude = 1 */
  int frame_ctr = 0;
#pragma omp parallel for
  for (int i = 0; i < FRAMES; i++) {
    char buff[100] = {0};
    sprintf(buff, "tt-%03d.png", i);
    FILE *fp = fopen(buff, "wb");
    assert(fp != NULL);
    printf("\r%s: %d/%d", buff, frame_ctr++, FRAMES - 1);
    fflush(stdout);

    png_structp png =
        png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    assert(png);
    png_infop info = png_create_info_struct(png);
    assert(info);
    int tem = setjmp(png_jmpbuf(png));
    assert(!tem);

    png_init_io(png, fp);
    png_set_IHDR(png, info, width, height, 8, PNG_COLOR_TYPE_RGBA,
                 PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
                 PNG_FILTER_TYPE_DEFAULT);
    png_write_info(png, info);

    png_bytep *row = NULL;
    row = malloc(sizeof(png_bytep) * height);
    assert(row != NULL);

    for (int y = 0; y < height; y++) {
      row[y] = malloc(sizeof(png_bytep) * width);
      assert(row[y] != NULL);
    }
    for (int y = 0; y < height; y++) {
      png_bytep line = row[y];
      for (int x = 0; x < width; x++) {
        double Px = x - (width / 2), Py = (height / 2) - y;
        struct vec3D ray = {Px, Py, Pz};
        double _magnitude = magnitude(ray);
        struct vec3D normalized_ray = {ray.x / _magnitude, ray.y / _magnitude,
                                       ray.z / _magnitude};
        /* RANDOM SPEED */
        struct object3D rotbox =
            rotate(box, (double)i * 3.0 / 180.0, M_PI * (double)i * 4.0 / 180.0,
                   (double)i * 7.0 / 180.0);
        struct object3D trbox = translate(rotbox, 0.0, 0.0, 100.0);
        free(rotbox.face);
        struct vec3D BmA;
        struct vec3D CmA;
        struct vec3D pvec;
        double illumination;
        double det, t = INFINITY; /* portable? */
        int intersect = 0;
        /* check if all triangle intersects */
        for (int triangle = 0; triangle < trbox.faces; triangle++) {
          BmA = (struct vec3D){.x = trbox.face[triangle].vertex[1].x -
                                    trbox.face[triangle].vertex[0].x,
                               .y = trbox.face[triangle].vertex[1].y -
                                    trbox.face[triangle].vertex[0].y,
                               .z = trbox.face[triangle].vertex[1].z -
                                    trbox.face[triangle].vertex[0].z};
          CmA = (struct vec3D){.x = trbox.face[triangle].vertex[2].x -
                                    trbox.face[triangle].vertex[0].x,
                               .y = trbox.face[triangle].vertex[2].y -
                                    trbox.face[triangle].vertex[0].y,
                               .z = trbox.face[triangle].vertex[2].z -
                                    trbox.face[triangle].vertex[0].z};
          pvec = (struct vec3D)cross_product(ray, CmA);
          double illum =
              dot_product(normalize(cross_product(BmA, CmA)), light_source);
          det = dot_product(BmA, pvec);

          double u = dot_product(pvec,
                                 (struct vec3D){
                                     .x = -trbox.face[triangle].vertex[0].x,
                                     .y = -trbox.face[triangle].vertex[0].y,
                                     .z = -trbox.face[triangle].vertex[0].z}) /
                     det;

          if (u < 0.0 || u > 1.0) {
            *(line + (x * 4) + 0) = 0xA0;
            *(line + (x * 4) + 1) = 0xA0;
            *(line + (x * 4) + 2) = 0xA0;
            *(line + (x * 4) + 3) = 0xFF;
            continue;
          }
          double v =
              dot_product(
                  cross_product(
                      (struct vec3D){.x = -trbox.face[triangle].vertex[0].x,
                                     .y = -trbox.face[triangle].vertex[0].y,
                                     .z = -trbox.face[triangle].vertex[0].z},
                      BmA),
                  ray) /
              det;
          if (v < 0.0 || u + v > 1.0) {
            *(line + (x * 4) + 0) = 0xA0;
            *(line + (x * 4) + 1) = 0xA0;
            *(line + (x * 4) + 2) = 0xA0;
            *(line + (x * 4) + 3) = 0xFF;
            continue;
          }
          double temp =
              dot_product(CmA, cross_product(
                                   (struct vec3D){
                                       .x = -trbox.face[triangle].vertex[0].x,
                                       .y = -trbox.face[triangle].vertex[0].y,
                                       .z = -trbox.face[triangle].vertex[0].z},
                                   BmA)) /
              det; /*
           double mm = magnitude((struct vec3D){normalized_ray.x*temp,
           normalized_ray.y*temp, normalized_ray.z*temp});*/
          if (t > temp) {
            illumination = illum;
            t = temp;
            intersect = 1;
          }
        }
        free(trbox.face);
        if (!intersect) {
          *(line + (x * 4) + 0) = 0xA0;
          *(line + (x * 4) + 1) = 0xA0;
          *(line + (x * 4) + 2) = 0xA0;
          *(line + (x * 4) + 3) = 0xFF;
          continue;
        }
        *(line + (x * 4) + 0) =
            (unsigned char)(((illumination + 1.0) / 2.0) * 255.0);
        *(line + (x * 4) + 1) =
            (unsigned char)(((illumination + 1.0) / 2.0) * 255.0);
        *(line + (x * 4) + 2) =
            (unsigned char)(((illumination + 1.0) / 2.0) * 255.0);
        *(line + (x * 4) + 3) = 0xFF;
      }
    }

    png_write_image(png, row);
    png_write_end(png, NULL);

    for (int y = 0; y < height; y++)
      free(row[y]);
    free(row);

    fclose(fp);
    png_destroy_write_struct(&png, &info);
  }
  putchar('\n');
  return 0;
}

faces.txt

[0].vertex = {[0].x = -20.0,
                    [0].y = -50.0,
                    [0].z = -20.0,
                    [1].x = -20.0,
                    [1].y = 50.0,
                    [1].z = -20.0,
                    [2].x = 20.0,
                    [2].y = 50.0,
                    [2].z = -20.0},
      [1].vertex = {[0].x = -20.0,
                    [0].y = -50.0,
                    [0].z = -20.0,
                    [1].x = 20.0,
                    [1].y = 50.0,
                    [1].z = -20.0,
                    [2].x = 20.0,
                    [2].y = -50.0,
                    [2].z = -20.0},
      [2].vertex = {[2].x = -20.0,
                    [2].y = -50.0,
                    [2].z = 20.0,
                    [1].x = -20.0,
                    [1].y = 50.0,
                    [1].z = 20.0,
                    [0].x = 20.0,
                    [0].y = 50.0,
                    [0].z = 20.0},
      [3].vertex = {[2].x = -20.0,
                    [2].y = -50.0,
                    [2].z = 20.0,
                    [1].x = 20.0,
                    [1].y = 50.0,
                    [1].z = 20.0,
                    [0].x = 20.0,
                    [0].y = -50.0,
                    [0].z = 20.0},
      [4].vertex = {[2].x = -20.0,
                    [2].y = -50.0,
                    [2].z = -20.0,
                    [1].x = -20.0,
                    [1].y = 50.0,
                    [1].z = -20.0,
                    [0].x = -20.0,
                    [0].y = -50.0,
                    [0].z = 20.0},
      [5].vertex = {[0].x = -20.0,
                    [0].y = -50.0,
                    [0].z = 20.0,
                    [1].x = -20.0,
                    [1].y = 50.0,
                    [1].z = 20.0,
                    [2].x = -20.0,
                    [2].y = 50.0,
                    [2].z = -20.0},
      [6].vertex = {[0].x = 20.0,
                    [0].y = -50.0,
                    [0].z = -20.0,
                    [1].x = 20.0,
                    [1].y = 50.0,
                    [1].z = -20.0,
                    [2].x = 20.0,
                    [2].y = -50.0,
                    [2].z = 20.0},
      [7].vertex = {[2].x = 20.0,
                    [2].y = -50.0,
                    [2].z = 20.0,
                    [1].x = 20.0,
                    [1].y = 50.0,
                    [1].z = 20.0,
                    [0].x = 20.0,
                    [0].y = 50.0,
                    [0].z = -20.0}