Ray casting with spheres.

.gitignore

@@ -27,3 +27,4 @@
 *.exe
 *.out
 *.app
+ray

Makefile

@@ -0,0 +1,20 @@
+TARGET = ray
+HEADERS = ray.hpp sphere.hpp figure.hpp
+OBJECTS = main.o sphere.o
+CXX = g++
+CXXFLAGS = -ansi -pedantic -Wall -g -DGLM_FORCE_RADIANS -fopenmp
+LDLIBS = -lm
+
+.PHONY: all
+all: $(TARGET)
+
+$(TARGET): $(OBJECTS) $(HEADERS)
+	$(CXX) $(CXXFLAGS) -o $@ $(OBJECTS) $(LDLIBS)
+
+main.o: main.cpp $(HEADERS)
+
+sphere.o: sphere.cpp $(HEADERS)
+
+.PHONY: clean
+clean:
+	$(RM) $(TARGET) *.o

figure.hpp

@@ -0,0 +1,26 @@
+#pragma once
+#ifndef FIGURE_HPP
+#define FIGURE_HPP
+
+#include <glm/vec3.hpp>
+
+#include "ray.hpp"
+
+class Figure {
+public:
+  vec3 color;
+
+  virtual ~Figure() { }
+
+  virtual bool intersect(Ray & r, float & t) const = 0;
+
+  virtual void set_color(float r, float g, float b) {
+    color = vec3(r, g, b);
+  }
+
+  virtual void set_color(vec3 rgb) {
+    color = vec3(rgb);
+  }
+};
+
+#endif

main.cpp

@@ -0,0 +1,176 @@
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+
+#include <omp.h>
+#include <glm/glm.hpp>
+
+#include "ray.hpp"
+#include "figure.hpp"
+#include "sphere.hpp"
+
+using namespace std;
+using namespace glm;
+
+static const vec3 BCKG_COLOR = vec3(0.16f, 0.66f, 0.72f);
+static const char * OUT_FILE = "output.ppm";
+
+static char * input_file;
+static int g_samples = 25;
+static float g_fov = 90.0f;
+static int g_w = 640;
+static int g_h = 480;
+static float g_aspect_ratio = static_cast<float>(g_w) / g_h;
+static vec3 ** image;
+
+float random01();
+vec2 sample_pixel(int i, int j);
+
+int main(int argc, char ** argv) {
+  FILE * out;
+  float t;
+  float _t;
+  Sphere * s;
+  Figure * _f;
+  Ray r;
+  vec2 sample;
+  vector<Figure *> figures;
+
+  if(argc < 2 || argc == 3 || argc > 6) {
+    cerr << "USAGE: " << argv[0] << " IN FILE [OUT FILE [HEIGHT WIDTH [SAMPLES [FIELD OF VIEW]]]]" << endl;
+    return EXIT_FAILURE;
+  }
+
+  input_file = argv[2];
+
+  if(argc >= 4) {
+    g_h = atoi(argv[2]);
+    if (g_h <= 0) {
+      cerr << "USAGE: " << argv[0] << " IN FILE [OUT FILE [HEIGHT WIDTH [SAMPLES [FIELD OF VIEW]]]]" << endl;
+      cerr << "HEIGHT must be positive" << endl;
+      return EXIT_FAILURE;
+    }
+
+    g_w = atoi(argv[3]);
+    if (g_w <= 0) {
+      cerr << "USAGE: " << argv[0] << " IN FILE [OUT FILE [HEIGHT WIDTH [SAMPLES [FIELD OF VIEW]]]]" << endl;
+      cerr << "WIDTH must be positive" << endl;
+      return EXIT_FAILURE;
+    }
+    
+    if(argc >= 5) {
+      g_samples = atoi(argv[4]);
+      if (g_samples <= 0) {
+	  cerr << "USAGE: " << argv[0] << " IN FILE [OUT FILE [HEIGHT WIDTH [SAMPLES [FIELD OF VIEW]]]]" << endl;
+	  cerr << "SAMPLES must be greater than 1" << endl;
+	  return EXIT_FAILURE;
+      }
+      
+      if(argc >= 6) {
+	g_fov = atof(argv[5]);
+	if (g_fov <= 0) {
+	  cerr << "USAGE: " << argv[0] << " IN FILE [OUT FILE [HEIGHT WIDTH [SAMPLES [FIELD OF VIEW]]]]" << endl;
+	  cerr << "FIELD OF VIEW must be greater than 1.0" << endl;
+	  return EXIT_FAILURE;
+	}
+      }
+    }
+  }
+
+  out = fopen(argc >= 2 ? argv[1] : OUT_FILE, "wb");
+
+  image = new vec3*[g_h];
+  for (int i = 0; i < g_h; i++) {
+    image[i] = new vec3[g_w];
+  }
+
+  s = new Sphere(1.0f, 1.0f, -2.0f, 0.5f);
+  s->set_color(1.0f, 0.0f, 0.0f);
+  figures.push_back(static_cast<Figure *>(s));
+
+  s = new Sphere(-1.0f, 1.0f, -2.0f, 0.5f);
+  s->set_color(0.0f, 1.0f, 0.0f);
+  figures.push_back(static_cast<Figure *>(s));
+
+  s = new Sphere(1.0f, -1.0f, -2.0f, 0.5f);
+  s->set_color(0.0f, 0.0f, 1.0f);
+  figures.push_back(static_cast<Figure *>(s));
+
+  s = new Sphere(-1.0f, -1.0f, -2.0f, 0.5f);
+  s->set_color(0.5f, 0.5f, 0.5f);
+  figures.push_back(static_cast<Figure *>(s));
+
+  s = new Sphere(0.0f, 0.0f, -2.0f, 1.0f);
+  s->set_color(1.0f, 1.0f, 0.0f);
+  figures.push_back(static_cast<Figure *>(s));
+
+#pragma omp parallel for schedule(dynamic, 1) private(r, sample, _f, t, _t)
+  for (int i = 0; i < g_h; i++) {
+    for (int j = 0; j < g_w; j++) {
+      for (int k = 0; k < g_samples; k++) {
+	sample = sample_pixel(i, j);
+	r = Ray(normalize(vec3(sample, -1.0f) - vec3(0.0f, 0.0f, 0.0f)), vec3(0.0f, 0.0f, 0.0f));
+	t = numeric_limits<float>::max();
+	_f = NULL;
+
+	for (size_t f = 0; f < figures.size(); f++) {
+	  if (figures[f]->intersect(r, _t) && _t < t) {
+	    t = _t;
+	    _f = figures[f];
+	  }
+	}
+
+	if (_f != NULL) {
+	  image[i][j] += vec3(_f->color);
+	} else {
+	  image[i][j] += vec3(BCKG_COLOR);
+	}
+      }
+
+      image[i][j] /= g_samples;
+    }
+  }
+
+  for (size_t f = 0; f < figures.size(); f++) {
+    delete static_cast<Sphere *>(figures[f]);
+  }
+  figures.clear();
+
+  fprintf(out, "P6 %d %d %d ", g_w, g_h, 255);
+  for (int i = 0; i < g_h; i++) {
+    for (int j = 0; j < g_w; j++) {
+      fputc(static_cast<int>(image[i][j].r * 255.0f), out);
+      fputc(static_cast<int>(image[i][j].g * 255.0f), out);
+      fputc(static_cast<int>(image[i][j].b * 255.0f), out);
+    }
+  }
+  fclose(out);
+
+  for (int i = 0; i < g_h; i++)
+    delete[] image[i];
+  delete[] image;
+
+  return EXIT_SUCCESS;
+}
+
+inline float random01() {
+  return static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
+}
+
+vec2 sample_pixel(int i, int j) {
+  float pxNDC;
+  float pyNDC;
+  float pxS;
+  float pyS;
+  pyNDC = (static_cast<float>(i) + random01()) / g_h;
+  pyS = (1.0f - (2.0f * pyNDC)) * tan(radians(g_fov) / 2);
+  pxNDC = (static_cast<float>(j) + random01()) / g_w;
+  pxS = (2.0f * pxNDC) - 1.0f;
+  pxS *= g_aspect_ratio * tan(radians(g_fov) / 2);
+
+  return vec2(pxS, pyS);
+}

ray.hpp

@@ -0,0 +1,19 @@
+#pragma once
+#ifndef RAY_HPP
+#define RAY_HPP
+
+#include <glm/vec3.hpp>
+
+using glm::vec3;
+
+class Ray {
+public:
+  vec3 m_direction;
+  vec3 m_origin;
+
+  Ray(): m_direction(vec3(0.0f, 0.0f, -1.0f)), m_origin(vec3(0.0f)) { }
+  Ray(float dx, float dy, float dz, float ox, float oy, float oz): m_direction(vec3(dx, dy, dz)), m_origin(vec3(ox, oy, oz)) {}
+  Ray(vec3 _d, vec3 _o): m_direction(_d), m_origin(_o) { }
+};
+
+#endif

sphere.cpp

@@ -0,0 +1,29 @@
+#include <cmath>
+
+#include "sphere.hpp"
+
+bool Sphere::intersect(Ray & r, float & t) const {
+  float d;
+
+  float a = (r.m_direction.x * r.m_direction.x) +
+    (r.m_direction.y * r.m_direction.y) +
+    (r.m_direction.z * r.m_direction.z);
+
+  float b = (2 * r.m_direction.x * (r.m_origin.x - m_center.x)) +
+    (2 * r.m_direction.y * (r.m_origin.y - m_center.y)) +
+    (2 * r.m_direction.z * (r.m_origin.z - m_center.z));
+
+  float c = (m_center.x * m_center.x) +
+    (m_center.y * m_center.y) +
+    (m_center.z * m_center.z) +
+    (r.m_origin.x * r.m_origin.x) +
+    (r.m_origin.y * r.m_origin.y) +
+    (r.m_origin.z * r.m_origin.z) -
+    2 * ((m_center.x * r.m_origin.x) + (m_center.y * r.m_origin.y) + (m_center.z * r.m_origin.z)) - (m_radius * m_radius);
+
+  d = (b * b) - (4 * a * c);
+
+  t = (-b - sqrt(d)) / (2 * a);
+
+  return d >= 0.0f;
+}

sphere.hpp

@@ -0,0 +1,27 @@
+#pragma once
+#ifndef SPHERE_HPP
+#define SPHERE_HPP
+
+#include <glm/vec3.hpp>
+
+#include "figure.hpp"
+
+using glm::vec3;
+
+class Sphere : public Figure {
+public:
+  vec3 m_center;
+  float m_radius;
+  
+  Sphere(): m_center(vec3(0.0f)), m_radius(0.5f) { }
+
+  Sphere(float x, float y, float z, float r): m_center(vec3(x, y, z)), m_radius(r) { }
+
+  Sphere(vec3 _c, float r): m_center(_c), m_radius(r) { }
+
+  virtual ~Sphere() { }
+
+  virtual bool intersect(Ray & r, float & t) const;
+};
+
+#endif