Minimal ray tracer for leaning purposes

ray.cc

#include <array>
#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>

struct Vec {
  float x, y, z;
  Vec(float vx, float vy, float vz) : x(vx), y(vy), z(vz) {}
  Vec operator+(Vec vec) { return {x + vec.x, y + vec.y, z + vec.z}; }
  Vec operator-(Vec vec) { return {x - vec.x, y - vec.y, z - vec.z}; }
  Vec operator*(float n) { return {x * n, y * n, z * n}; }
  Vec unit() { return Vec(x, y, z) * (1 / this->length()); }
  // dot product
  float operator%(Vec vec) { return x * vec.x + y * vec.y + z * vec.z; }
  float length() { return sqrtf(*this % *this); }
};

struct Sphere {
  Vec center;
  float color;
  float radius;
  float intersect(Vec origin, Vec direction) {
    Vec p = origin - this->center;
    float a = direction % direction;
    float b = (p % direction) * 2;
    float c = (p % p) - (this->radius * this->radius);
    float d = b * b - 4 * a * c;
    if (d < 0) {
      return NAN;
    }
    float sqd = sqrtf(d);
    float distance = (-b - sqd) / (2.f * a);
    if (distance > .1f) {
      return distance;
    }
    distance = (-b + sqd) / (2.f * a);
    if (distance > .1f) {
      return distance;
    }
    return NAN;
  }
};

struct World {
  std::vector<Sphere> spheres;
  std::vector<Sphere> lights;
};

float trace(World world, Vec origin, Vec direction) {
  int index = -1;
  float distance = NAN;
  for (int i = 0; i < world.spheres.size(); ++i) {
    float d = world.spheres[i].intersect(origin, direction);
    if (!std::isnan(d) && (index < 0 || d < distance)) {
      distance = d;
      index = i;
    }
  }
  if (index < 0) {
    return 1.f - direction.y;
  }

  Vec p = origin + direction * distance;
  Vec n = (p - world.spheres[index].center).unit();
  float c = world.spheres[index].color * .1f;
  for (auto light : world.lights) {
    Vec l = (light.center - p).unit();
    int shadow = 0;
    for (auto sphere : world.spheres) {
      if (!std::isnan(sphere.intersect(p, l))) {
        shadow = 1;
      }
    }
    if (!shadow) {
      float df = std::max(0.f, (l % n) * 0.7f);
      float sp = powf(fmax(0.f, (l % n)), 70.f) * 0.4f;
      c = c + world.spheres[index].color * light.color * df + sp;
    }
  }
  return c;
}

void render_pgm_stereo(World world, std::string filename, int width,
                       int height) {
  std::ofstream f(filename);
  f << "P2" << std::endl << (width * 2) << " " << height << " 255" << std::endl;
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      float c = trace(world, {0, 1, 5},
                      Vec(x - width / 2, height / 2 - y, -height).unit());
      f << ((int)(c * 255)) << " ";
    }
    for (int x = 0; x < width; x++) {
      float c = trace(world, {0.5, 1, 5},
                      Vec(x - width / 2, height / 2 - y, -height).unit());
      f << ((int)(c * 255)) << " ";
    }
  }
}

void render_pgm(World world, std::string filename, int width, int height) {
  std::ofstream f(filename);
  f << "P2" << std::endl << width << " " << height << " 255" << std::endl;
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      float c = trace(world, {0, 1, 5},
                      Vec(x - width / 2, height / 2 - y, -height).unit());
      f << ((int)(c * 255)) << " ";
    }
  }
}

void render_tty(World world, int width, int height) {
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      float c = trace(world, {0, 1, 5},
                      Vec(x - width / 2, height / 2 - y, -height).unit());
      char pixel = " .:-=+*#%@$"[std::max(std::min((int)(c * 10), 10), 0)];
      std::cout << pixel << pixel;
    }
    std::cout << std::endl;
  }
}

int main() {
  World world = {
      // spheres
      {
          {{0, -1000, 0}, 0.001, 1000},
          {{-2, 1, -2}, 1, 1},
          {{0, 1, 0}, 0.5, 1},
          {{2, 1, -1}, 0.1, 1},
      },
      // lights
      {
          {{0, 100, 0}, .4, 0},
          {{100, 100, 200}, .5, 0},
          {{-100, 300, 100}, .1, 0},
      },
  };

  render_pgm(world, "ray_mono.pgm", 600, 600);
  render_pgm_stereo(world, "ray.pgm", 300, 200);
  render_tty(world, 40, 25);
  return 0;
}

Hey @zserge, nice blog post!

I think you have the definition of the cross product wrong. Given vectors A(1, 0, 0) and B(0, 1, 0), the cross product should be C(0, 0, 1), but your definition return {x * vec.x, y * vec.y, z * vec.z} gives C(0, 0, 0). The correct definition is more like this:

(source: https://en.wikipedia.org/wiki/Cross_product)

Also, the cross product is entirelly unnecessary for your ray tracer (which is why it doesn't matter that its definiton is incorrect), as can be easily verified by commenting out the definition of Vec operator*(Vec vec).