3D scene

main: main.c buf.c
	$(CC) -o $@ $^ -lm

main: buf.h

main.pbm: main
	./main

.PHONY: clean

clean:
	$(RM) main main.pbm

#include "buf.h"
#include <stdio.h>
#include <stdlib.h>

void
file_out(const struct buf *b, const char *filename)
{
   const int width  = b->width;
   const int height = b->height;
   const unsigned char *buffer = b->bytes;

   FILE *f = fopen(filename, "w");

   fprintf(f, "P1\n%d %d", width, height);
   for (int i = 0; i < width * height; i++) {
      fputc(i % 32 == 0 ? '\n' : ' ', f);
      fputc(buffer[i] ? '1' : '0', f); }
   fputc('\n', f);

   fclose(f);
}

static void swap(int *a, int *b) { const int t = *a; *a = *b, *b = t; }

void
line(const struct buf *b, int x0, int y0, int x1, int y1)
{
   if (y0 == y1) {
      if (x0 > x1) swap(&x0, &x1);
      for (int x = x0; x <= x1; x++)
         pset(b, x, y0, 1); }
   else if (x0 == x1) {
      if (y0 > y1) swap(&y0, &y1);
      for (int y = y0; y <= y1; y++)
         pset(b, x0, y, 1); }
   else if (abs(x1 - x0) > abs(y1 - y0)) {
      if (x0 > x1) swap(&x0, &x1), swap(&y0, &y1);
      for (int x = x0; x <= x1; x++) {
         const int y = y0 + (double)(x - x0) / (x1 - x0) * (y1 - y0);
         pset(b, x, y, 1); } }
   else {
      if (y0 > y1) swap(&x0, &x1), swap(&y0, &y1);
      for (int y = y0; y <= y1; y++) {
         const int x = x0 + (double)(y - y0) / (y1 - y0) * (x1 - x0);
         pset(b, x, y, 1); } }
}

struct buf {
   int width;
   int height;
   unsigned char *bytes;
};

void file_out(const struct buf *, const char *);
void line(const struct buf *b, int, int, int, int);

static void
pset(const struct buf *b, const int x, const int y, const int val)
{
   const int i = x + y * b->width;

   if (i < 0 || i > b->width * b->height) return;
   if (x < 0 || x >= b->width)  return;
   if (y < 0 || y >= b->height) return;

   b->bytes[i] = val ? 1 : 0;
}

#include "buf.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h> // XXX Debug
#define W 640
#define H 480
#define DPDU 600 // dots per distance unit (DU)

// config
static const double camera_x =  0.3;
static const double camera_y = -0.2;
static const double camera_z = -3.0;

typedef struct { double x, y, z ; } vec;
typedef struct { double x, y    ; } vec2d;

// distance from eye to screen is 1 DU.
double screen_x(double x_3d, double z_3d) { return x_3d / z_3d; }
double screen_y(double y_3d, double z_3d) { return y_3d / z_3d; }

int
main(void)
{
   unsigned char bytes[W * H] = { 0 };
   const struct buf b = {
      .width  = W,
      .height = H,
      .bytes  = bytes,
   };

#if 0
   const vec dots[] = {
      // distance = 1.0 (same as the screen)
      { -0.9,  0.9,  1.0 }, {  0.0,  0.9,  1.0 }, {  0.9,  0.9,  1.0 },
      { -0.9,  0.0,  1.0 }, {  0.0,  0.0,  1.0 }, {  0.9,  0.0,  1.0 },
      { -0.9, -0.9,  1.0 }, {  0.0, -0.9,  1.0 }, {  0.9, -0.9,  1.0 },
      // distance = 2.0
      { -1.0,  1.0,  2.0 }, {  0.0,  1.0,  2.0 }, {  1.0,  1.0,  2.0 },
      { -1.0,  0.0,  2.0 }, {  0.0,  0.0,  2.0 }, {  1.0,  0.0,  2.0 },
      { -1.0, -1.0,  2.0 }, {  0.0, -1.0,  2.0 }, {  1.0, -1.0,  2.0 },
      // distance = 3.0
      { -1.0,  1.0,  3.0 }, {  0.0,  1.0,  3.0 }, {  1.0,  1.0,  3.0 },
      { -1.0,  0.0,  3.0 }, {  0.0,  0.0,  3.0 }, {  1.0,  0.0,  3.0 },
      { -1.0, -1.0,  3.0 }, {  0.0, -1.0,  3.0 }, {  1.0, -1.0,  3.0 },
      // distance = 4.0
      { -1.0,  1.0,  4.0 }, {  0.0,  1.0,  4.0 }, {  1.0,  1.0,  4.0 },
      { -1.0,  0.0,  4.0 }, {  0.0,  0.0,  4.0 }, {  1.0,  0.0,  4.0 },
      { -1.0, -1.0,  4.0 }, {  0.0, -1.0,  4.0 }, {  1.0, -1.0,  4.0 },
   };
   const size_t n_dots = sizeof dots / sizeof *dots;

   for (int i = 0; i < n_dots; i++) {
      const vec *v = &dots[i];
      const double sx = screen_x(v->x, v->z);
      const double sy = screen_y(v->y, v->z);
      const int bx = sx * DPDU + W / 2;
      const int by = sy * DPDU + H / 2;

      printf("(%3d, %3d) <- (%5.2f, %5.2f)\n", bx, by, sx, sy); // XXX Debug

#define LINE(dx0, dy0, dx1, dy1) \
      line(&b, bx + dx0, H - (by + dy0), bx + dx1, H - (by + dy1))
      LINE(-4, -4, +4, +4); // "\"
      LINE(+4, -4, -4, +4); // "/"
      LINE( 0, -4,  0, +4); // "|"
      LINE(-4,  0, +4,  0); // "-"

      pset(&b, bx, H - by, 1); }
#endif

   const vec2d lines[] = {
      { -1.0, -1.0 },
      {  1.0, -1.0 },
      {  1.0,  1.0 },
      { -1.0,  1.0 },
      { -1.0, -1.0 },
   };

   for (int i = -2; i < 16; i++) {
      const double z = 1.0 * (i + 1) - camera_z;

      for (int j = 0; j < 4; j++) {
         const vec2d f = lines[j]; // from
         const vec2d t = lines[j + 1]; // to
         const double sx0 = screen_x(f.x - camera_x, z);
         const double sy0 = screen_y(f.y - camera_y, z);
         const double sx1 = screen_x(t.x - camera_x, z);
         const double sy1 = screen_y(t.y - camera_y, z);
         const int bx0 = sx0 * DPDU + W / 2;
         const int by0 = sy0 * DPDU + H / 2;
         const int bx1 = sx1 * DPDU + W / 2;
         const int by1 = sy1 * DPDU + H / 2;

         line(&b, bx0, H - by0, bx1, H - by1); } }

#if 1
   const vec2d floor_lines[] = { // use y as z
      { -0.6, -0.6 },
      {  0.6, -0.6 },
      {  0.6,  0.6 },
      { -0.6,  0.6 },
      { -0.6, -0.6 },
   };

   for (int i = 0; i < 4; i++) {
      const vec2d f = floor_lines[i];
      const vec2d t = floor_lines[i + 1];
      const double z0 = f.y - camera_z;
      const double z1 = t.y - camera_z;
      const double sx0 = screen_x(f.x  - camera_x, z0);
      const double sy0 = screen_y(-1.0 - camera_y, z0);
      const double sx1 = screen_x(t.x  - camera_x, z1);
      const double sy1 = screen_y(-1.0 - camera_y, z1);
      const int bx0 = sx0 * DPDU + W / 2;
      const int by0 = sy0 * DPDU + H / 2;
      const int bx1 = sx1 * DPDU + W / 2;
      const int by1 = sy1 * DPDU + H / 2;

      line(&b, bx0, H - by0, bx1, H - by1); }
#endif

#if 1
   for (int i = 0; i < 24; i++) {
      const double pi = 3.14;
      const double th0 = 2 * pi * (i + 0) / 24.0;
      const double th1 = 2 * pi * (i + 4) / 24.0;
      const vec2d f = { .x = 0.6 * cos(th0), .y = 0.6 * sin(th0), };
      const vec2d t = { .x = 0.6 * cos(th1), .y = 0.6 * sin(th1), };
      const double z0 = f.y - camera_z;
      const double z1 = t.y - camera_z;
      const double sx0 = screen_x(f.x  - camera_x, z0);
      const double sy0 = screen_y(-1.0 - camera_y, z0);
      const double sx1 = screen_x(t.x  - camera_x, z1);
      const double sy1 = screen_y(-1.0 - camera_y, z1);
      const int bx0 = sx0 * DPDU + W / 2;
      const int by0 = sy0 * DPDU + H / 2;
      const int bx1 = sx1 * DPDU + W / 2;
      const int by1 = sy1 * DPDU + H / 2;

      line(&b, bx0, H - by0, bx1, H - by1); }
#endif

#if 1
   for (int h = 0; h < 4; h++) {
      const double r = 1.0 + h * 0.8; // radius

      const int N = 24;
      for (int i = 0; i < N; i++) {
         const double pi = 3.14;
         const double th0 = 2 * pi * (i + 0) / N;
         const double th1 = 2 * pi * (i + 1) / N;
         const vec2d f = { .x = r * cos(th0), .y = r * sin(th0), };
         const vec2d t = { .x = r * cos(th1), .y = r * sin(th1), };
         const double z0 = f.y - camera_z;
         const double z1 = t.y - camera_z;
         const double sx0 = screen_x(f.x  - camera_x, z0);
         const double sy0 = screen_y(-1.0 - camera_y, z0);
         const double sx1 = screen_x(t.x  - camera_x, z1);
         const double sy1 = screen_y(-1.0 - camera_y, z1);
         const int bx0 = sx0 * DPDU + W / 2;
         const int by0 = sy0 * DPDU + H / 2;
         const int bx1 = sx1 * DPDU + W / 2;
         const int by1 = sy1 * DPDU + H / 2;
         // const int bx0 = sx0 / 10 * DPDU + W / 2; // XXX Debug
         // const int by0 = sy0 / 10 * DPDU + H / 2;
         // const int bx1 = sx1 / 10 * DPDU + W / 2;
         // const int by1 = sy1 / 10 * DPDU + H / 2;

         if (z0 < 1.0 || z1 < 1.0) continue;

         // XXX Debug
         const int bxd = bx0 - bx1, byd = by0 - by1;
         const double d = sqrt((double)bxd * bxd + (double)byd * byd);
         if (d > 1000.0) {
            printf("(%d %d)-(%d %d) [%d %d %.0f]:[%.2f %.2f]\n",
                  bx0, by0, bx1, by1, h, i, d, z0, z1);
            continue; }

         line(&b, bx0, H - by0, bx1, H - by1); }
   }
#endif

#if 1
   {
      // hidden line removal - painters algorithm
      int xmax[H], xmin[H];
      for (int i = 0; i < H; i++) xmax[i] = -W, xmin[i] = W;

      const int x_centre = screen_x(-camera_x, -camera_z) * DPDU;
      const double pi = 3.14;
      const int N = 80;

      for (int h = 0; h < 8 * N; h++) {
         const double h0 = (double)h / N;
         const double y0 = -1.0 + h0 * 0.30;
         const double r0 =  0.6 - h0 * 0.07; // radius
         const double phase0 =    h0 * 10.0 / 360.0 * 2 * pi; // rad

         for (int i = 0; i < 3; i++) {
            const double th0 = 2 * pi * i / 3.0 + phase0;
            const vec2d f = { .x = r0 * cos(th0), .y = r0 * sin(th0), };
            const double z0 = f.y - camera_z;
            const double sx0 = screen_x(f.x - camera_x, z0);
            const double sy0 = screen_y(y0  - camera_y, z0);
            const int bx0 = sx0 * DPDU - x_centre;
            const int by0 = sy0 * DPDU + H / 2;

            if (by0 >= H || by0 < 0) continue;
            if (xmax[by0] < bx0) xmax[by0] = bx0;
            if (xmin[by0] > bx0) xmin[by0] = bx0; } }

      for (int i = 0; i < H; i++)
         for (int j = xmin[i]; j < xmax[i]; j++)
            pset(&b, W / 2 + x_centre + j, H - i, 0);

#if 0 // XXX Debug
      for (int i = 0; i < H; i++)
         printf("%6d  [%d %d]\n", i, xmin[i], xmax[i]);
#endif
   }

   for (int h = 0; h < 8; h++) {
      const double y = -1.0 + h * 0.30;
      const double r =  0.6 - h * 0.07; // radius
      const double phase =    h * 10.0 / 360.0 * 2 * 3.14; // rad

      for (int i = 0; i < 3; i++) {
         const double pi = 3.14;
         const double th0 = 2 * pi * (i + 0) / 3.0 + phase;
         const double th1 = 2 * pi * (i + 1) / 3.0 + phase;
         const vec2d f = { .x = r * cos(th0), .y = r * sin(th0), };
         const vec2d t = { .x = r * cos(th1), .y = r * sin(th1), };
         const double z0 = f.y - camera_z;
         const double z1 = t.y - camera_z;
         const double sx0 = screen_x(f.x  - camera_x, z0);
         const double sy0 = screen_y(y    - camera_y, z0);
         const double sx1 = screen_x(t.x  - camera_x, z1);
         const double sy1 = screen_y(y    - camera_y, z1);
         const int bx0 = sx0 * DPDU + W / 2;
         const int by0 = sy0 * DPDU + H / 2;
         const int bx1 = sx1 * DPDU + W / 2;
         const int by1 = sy1 * DPDU + H / 2;

         line(&b, bx0, H - by0, bx1, H - by1); }
   }

   for (int h = 0; h < 8; h++) {
      const int h0 = h, h1 = h + 1;
      const double pi = 3.14;
      const double y0 = -1.0 + h0 * 0.30;
      const double y1 = -1.0 + h1 * 0.30;
      const double r0 =  0.6 - h0 * 0.07; // radius
      const double r1 =  0.6 - h1 * 0.07; // radius
      const double phase0 =   h0 * 10.0 / 360.0 * 2 * pi; // rad
      const double phase1 =   h1 * 10.0 / 360.0 * 2 * pi; // rad

      for (int i = 0; i < 3; i++) {
         const double th0 = 2 * pi * i / 3.0 + phase0;
         const double th1 = 2 * pi * i / 3.0 + phase1;
         const vec2d f = { .x = r0 * cos(th0), .y = r0 * sin(th0), };
         const vec2d t = { .x = r1 * cos(th1), .y = r1 * sin(th1), };
         const double z0 = f.y - camera_z;
         const double z1 = t.y - camera_z;
         const double sx0 = screen_x(f.x  - camera_x, z0);
         const double sy0 = screen_y(y0   - camera_y, z0);
         const double sx1 = screen_x(t.x  - camera_x, z1);
         const double sy1 = screen_y(y1   - camera_y, z1);
         const int bx0 = sx0 * DPDU + W / 2;
         const int by0 = sy0 * DPDU + H / 2;
         const int bx1 = sx1 * DPDU + W / 2;
         const int by1 = sy1 * DPDU + H / 2;

         line(&b, bx0, H - by0, bx1, H - by1); }
   }
#endif

   file_out(&b, "main.pbm");
   return 0;
}