triangle.c

#include <assert.h>
#include <math.h>
#include <png.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define FRAMES 300
struct vec3D {
  double x, y, z;
};

struct triangle3D {
  struct vec3D vertex[3];
};

struct object3D {
  int faces;
  struct triangle3D *face;
};

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, c, r;
    Pz = 150.0;
  width = 1000, height = 720;
  struct triangle3D face = {.vertex = {[0].x = -20.0,
                                       [0].y = 0.0,
                                       [0].z = 0.0,
                                       [1].x = 0.0,
                                       [1].y = 90.0,
                                       [1].z = 0.0,
                                       [2].x = 170.0,
                                       [2].y = -200.0,
                                       [2].z = 0.0}};
  struct object3D box = {.faces = 1, .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);
    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, M_PI*(double)i*3.0/180.0, 2+(double)i*4/181.0, (double)i/180.0);
        struct object3D trbox = translate(rotbox, 0.0, 0.0, 100.0);
        free(rotbox.face);
        /* compute surface normal */
        struct vec3D BmA = {
            .x = trbox.face[0].vertex[1].x - trbox.face[0].vertex[0].x,
            .y = trbox.face[0].vertex[1].y - trbox.face[0].vertex[0].y,
            .z = trbox.face[0].vertex[1].z - trbox.face[0].vertex[0].z};
        struct vec3D CmA = {
            .x = trbox.face[0].vertex[2].x - trbox.face[0].vertex[0].x,
            .y = trbox.face[0].vertex[2].y - trbox.face[0].vertex[0].y,
            .z = trbox.face[0].vertex[2].z - trbox.face[0].vertex[0].z};
        struct vec3D pvec = cross_product(ray, CmA);
        double illumination =
            dot_product(normalize(cross_product(BmA, CmA)), light_source);
        double det = dot_product(BmA, pvec);

        double u =
            dot_product(pvec, (struct vec3D){.x = -trbox.face[0].vertex[0].x,
                                             .y = -trbox.face[0].vertex[0].y,
                                             .z = -trbox.face[0].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;
          free(trbox.face);
          continue;
        }

        double v =
            dot_product(
                cross_product((struct vec3D){.x = -trbox.face[0].vertex[0].x,
                                             .y = -trbox.face[0].vertex[0].y,
                                             .z = -trbox.face[0].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;
          free(trbox.face);
          continue;
        }
        free(trbox.face);
        *(line + (x * 4) + 0) = (unsigned char)(fabs(illumination) * 255.0);
        *(line + (x * 4) + 1) = (unsigned char)(fabs(illumination) * 255.0);
        *(line + (x * 4) + 2) = (unsigned char)(fabs(illumination) * 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;
}