depth map raytracing

.gitpod.yml

@@ -0,0 +1,12 @@
+image:
+  file: Dockerfile
+tasks:
+- command: >
+    mkdir --parents build &&
+    cd build &&
+    cmake .. &&
+    make &&
+    ./tinyraycaster &&
+    pnmtopng out.ppm > out.png &&
+    open out.png &&
+    cd ..

Dockerfile

@@ -0,0 +1,6 @@
+FROM gitpod/workspace-full
+
+USER root
+# add your tools here
+RUN apt-get update && apt-get install -y \
+  netpbm

geometry.h

@@ -3,83 +3,41 @@
 #include <cmath>
 #include <vector>
 #include <cassert>
-#include <iostream>
 
-template <size_t DIM, typename T> struct vec {
+template <typename T> struct vec3 {
-    vec() { for (size_t i=DIM; i--; data_[i] = T()); }
-          T& operator[](const size_t i)       { assert(i<DIM); return data_[i]; }
-    const T& operator[](const size_t i) const { assert(i<DIM); return data_[i]; }
-private:
-    T data_[DIM];
-};
-
-typedef vec<2, float> Vec2f;
-typedef vec<3, float> Vec3f;
-typedef vec<3, int  > Vec3i;
-typedef vec<4, float> Vec4f;
-
-template <typename T> struct vec<2,T> {
-    vec() : x(T()), y(T()) {}
-    vec(T X, T Y) : x(X), y(Y) {}
-    template <class U> vec<2,T>(const vec<2,U> &v);
-          T& operator[](const size_t i)       { assert(i<2); return i<=0 ? x : y; }
-    const T& operator[](const size_t i) const { assert(i<2); return i<=0 ? x : y; }
-    T x,y;
-};
-
-template <typename T> struct vec<3,T> {
-    vec() : x(T()), y(T()), z(T()) {}
-    vec(T X, T Y, T Z) : x(X), y(Y), z(Z) {}
           T& operator[](const size_t i)       { assert(i<3); return i<=0 ? x : (1==i ? y : z); }
     const T& operator[](const size_t i) const { assert(i<3); return i<=0 ? x : (1==i ? y : z); }
-    float norm() { return std::sqrt(x*x+y*y+z*z); }
+    vec3<T> & normalize(T l=1) { *this = (*this)*(l/std::sqrt((*this)*(*this))); return *this; }
-    vec<3,T> & normalize(T l=1) { *this = (*this)*(l/norm()); return *this; }
     T x,y,z;
 };
 
-template <typename T> struct vec<4,T> {
+template<typename T> T operator*(const vec3<T>& lhs, const vec3<T>& rhs) {
-    vec() : x(T()), y(T()), z(T()), w(T()) {}
-    vec(T X, T Y, T Z, T W) : x(X), y(Y), z(Z), w(W) {}
-          T& operator[](const size_t i)       { assert(i<4); return i<=0 ? x : (1==i ? y : (2==i ? z : w)); }
-    const T& operator[](const size_t i) const { assert(i<4); return i<=0 ? x : (1==i ? y : (2==i ? z : w)); }
-    T x,y,z,w;
-};
-
-template<size_t DIM,typename T> T operator*(const vec<DIM,T>& lhs, const vec<DIM,T>& rhs) {
     T ret = T();
-    for (size_t i=DIM; i--; ret+=lhs[i]*rhs[i]);
+    for (size_t i=3; i--; ret+=lhs[i]*rhs[i]);
     return ret;
 }
 
-template<size_t DIM,typename T>vec<DIM,T> operator+(vec<DIM,T> lhs, const vec<DIM,T>& rhs) {
+template<typename T> vec3<T> operator+(vec3<T> lhs, const vec3<T>& rhs) {
-    for (size_t i=DIM; i--; lhs[i]+=rhs[i]);
+    for (size_t i=3; i--; lhs[i]+=rhs[i]);
     return lhs;
 }
 
-template<size_t DIM,typename T>vec<DIM,T> operator-(vec<DIM,T> lhs, const vec<DIM,T>& rhs) {
+template<typename T> vec3<T> operator-(vec3<T> lhs, const vec3<T>& rhs) {
-    for (size_t i=DIM; i--; lhs[i]-=rhs[i]);
+    for (size_t i=3; i--; lhs[i]-=rhs[i]);
     return lhs;
 }
 
-template<size_t DIM,typename T,typename U> vec<DIM,T> operator*(const vec<DIM,T> &lhs, const U& rhs) {
+template<typename T,typename U> vec3<T> operator*(const vec3<T> &lhs, const U& rhs) {
-    vec<DIM,T> ret;
+    vec3<T> ret;
-    for (size_t i=DIM; i--; ret[i]=lhs[i]*rhs);
+    for (size_t i=3; i--; ret[i]=lhs[i]*rhs);
     return ret;
 }
 
-template<size_t DIM,typename T> vec<DIM,T> operator-(const vec<DIM,T> &lhs) {
+template<typename T> vec3<T> operator-(const vec3<T> &lhs) {
     return lhs*T(-1);
 }
 
-template <typename T> vec<3,T> cross(vec<3,T> v1, vec<3,T> v2) {
+typedef vec3<float> Vec3f;
-    return vec<3,T>(v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z, v1.x*v2.y - v1.y*v2.x);
-}
 
-template <size_t DIM, typename T> std::ostream& operator<<(std::ostream& out, const vec<DIM,T>& v) {
-    for(unsigned int i=0; i<DIM; i++) {
-        out << v[i] << " " ;
-    }
-    return out ;
-}
 #endif //__GEOMETRY_H__
 

tinyraytracer.cpp

@@ -1,100 +1,94 @@
-#include <limits>
+#define _USE_MATH_DEFINES
 #include <cmath>
+#include <limits>
 #include <iostream>
 #include <fstream>
 #include <vector>
+#include <numeric>
 #include "geometry.h"
 
-struct Material {
-    Material(const Vec3f &color) : diffuse_color(color) {}
-    Material() : diffuse_color() {}
-    Vec3f diffuse_color;
-};
-
 struct Sphere {
     Vec3f center;
     float radius;
-    Material material;
-
-    Sphere(const Vec3f &c, const float &r, const Material &m) : center(c), radius(r), material(m) {}
-
-    bool ray_intersect(const Vec3f &orig, const Vec3f &dir, float &t0) const {
-        Vec3f L = center - orig;
-        float tca = L*dir;
-        float d2 = L*L - tca*tca;
-        if (d2 > radius*radius) return false;
-        float thc = sqrtf(radius*radius - d2);
-        t0       = tca - thc;
-        float t1 = tca + thc;
-        if (t0 < 0) t0 = t1;
-        if (t0 < 0) return false;
-        return true;
-    }
 };
 
-bool scene_intersect(const Vec3f &orig, const Vec3f &dir, const std::vector<Sphere> &spheres, Vec3f &hit, Vec3f &N, Material &material) {
+bool sphere_ray_intersect(const Sphere &sphere, const Vec3f &orig, const Vec3f &dir, float &t0) {
+    Vec3f L = sphere.center - orig;
+    float tca = L*dir;
+    float d2 = L*L - tca*tca;
+    if (d2 > pow(sphere.radius, 2)) return false;
+    float thc = sqrtf(pow(sphere.radius, 2) - d2);
+    t0       = tca - thc;
+    float t1 = tca + thc;
+    if (t0 < 0) t0 = t1;
+    if (t0 < 0) return false;
+    return true;
+}
+
+float scene_intersect(const Vec3f &orig, const Vec3f &dir, const std::vector<Sphere> &spheres) {
     float spheres_dist = std::numeric_limits<float>::max();
-    for (size_t i=0; i < spheres.size(); i++) {
+    for (size_t i=0; i<spheres.size(); i++) {
         float dist_i;
-        if (spheres[i].ray_intersect(orig, dir, dist_i) && dist_i < spheres_dist) {
+        if (sphere_ray_intersect(spheres[i], orig, dir, dist_i) && dist_i < spheres_dist) {
             spheres_dist = dist_i;
-            hit = orig + dir*dist_i;
-            N = (hit - spheres[i].center).normalize();
-            material = spheres[i].material;
         }
     }
-    return spheres_dist<1000;
+
+    float checkerboard_dist = std::numeric_limits<float>::max();
+    if (fabs(dir.y)>1e-3) {
+        float d = -(orig.y+4)/dir.y; // the checkerboard plane has equation y = -4
+        Vec3f pt = orig + dir*d;
+        if (d>0 && fabs(pt.x)<10 && pt.z<-10 && pt.z>-30) {
+            checkerboard_dist = d;
+        }
+    }
+
+    return std::min(spheres_dist, checkerboard_dist);
 }
 
-Vec3f cast_ray(const Vec3f &orig, const Vec3f &dir, const std::vector<Sphere> &spheres) {
+void compute_depthmap(const size_t width, const size_t height, const float fov, const float far, const std::vector<Sphere> &spheres, std::vector<float> &zbuffer) {
-    Vec3f point, N;
+#pragma omp parallel for
-    Material material;
+    for (size_t j=0; j<height; j++) {
-
+        for (size_t i=0; i<width; i++) {
-    if (!scene_intersect(orig, dir, spheres, point, N, material)) {
+            float dir_x =  (i + 0.5) - ( width/2.);
-        return Vec3f(0.2, 0.7, 0.8); // background color
+            float dir_y = -(j + 0.5) + (height/2.);   // this flips the image at the same time
-    }
+            float dir_z = -(height/(2.*tan(fov/2.)));
-
+            Vec3f dir = Vec3f{dir_x, dir_y, dir_z}.normalize();
-    return material.diffuse_color;
+            zbuffer[i+j*width] = scene_intersect(Vec3f{0,0,0}, dir, spheres);
-}
-
-void render(const std::vector<Sphere> &spheres) {
-    const int width    = 1024;
-    const int height   = 768;
-    const int fov      = M_PI/2.;
-    std::vector<Vec3f> framebuffer(width*height);
-
-    #pragma omp parallel for
-    for (size_t j = 0; j<height; j++) {
-        for (size_t i = 0; i<width; i++) {
-            float x =  (2*(i + 0.5)/(float)width  - 1)*tan(fov/2.)*width/(float)height;
-            float y = -(2*(j + 0.5)/(float)height - 1)*tan(fov/2.);
-            Vec3f dir = Vec3f(x, y, -1).normalize();
-            framebuffer[i+j*width] = cast_ray(Vec3f(0,0,0), dir, spheres);
         }
     }
 
-    std::ofstream ofs; // save the framebuffer to file
+    float minval =  std::numeric_limits<float>::max(); // clamp the zbuffer by the far plane and normalize it between 0 and 1
-    ofs.open("./out.ppm");
+    float maxval = -std::numeric_limits<float>::max();
-    ofs << "P6\n" << width << " " << height << "\n255\n";
+    for (size_t i=0; i<height*width; ++i) {
-    for (size_t i = 0; i < height*width; ++i) {
+        minval = std::min(minval, zbuffer[i]);
-        for (size_t j = 0; j<3; j++) {
+        maxval = std::max(maxval, std::min(zbuffer[i], far));
-            ofs << (char)(255 * std::max(0.f, std::min(1.f, framebuffer[i][j])));
-        }
     }
-    ofs.close();
+    for (size_t i=0; i<height*width; ++i)
+        zbuffer[i] = 1-(std::min(zbuffer[i],far)-minval)/(maxval-minval);
 }
 
 int main() {
-    Material      ivory(Vec3f(0.4, 0.4, 0.3));
+    std::vector<Sphere> spheres = { {{-3,0,-16}, 2}, {{-1,-1.5,-12}, 2}, {{1.5,-.5,-18}, 3}, {{7,5,-18}, 4} };
-    Material red_rubber(Vec3f(0.3, 0.1, 0.1));
+    const size_t width  = 1024;
+    const size_t height = 768;
+    const float fov  = M_PI/3.;
+    std::vector<float> zbuffer(width*height);
+    compute_depthmap(width, height, fov, 23.f, spheres, zbuffer);
 
-    std::vector<Sphere> spheres;
+    std::vector<unsigned char> framebuffer(width*height*3);
-    spheres.push_back(Sphere(Vec3f(-3,    0,   -16), 2,      ivory));
+    for (size_t j=0; j<height; j++) {
-    spheres.push_back(Sphere(Vec3f(-1.0, -1.5, -12), 2, red_rubber));
+        for (size_t i=0; i<width; i++) {
-    spheres.push_back(Sphere(Vec3f( 1.5, -0.5, -18), 3, red_rubber));
+            for (size_t d=0; d<3; d++) {
-    spheres.push_back(Sphere(Vec3f( 7,    5,   -18), 4,      ivory));
+                framebuffer[(i+j*width)*3 + d] = 255*zbuffer[i+j*width];
+            }
+        }
+    }
 
-    render(spheres);
+    std::ofstream ofs("./out.ppm"); // save the framebuffer to file
+    ofs << "P6\n" << width << " " << height << "\n255\n";
+    ofs.write(reinterpret_cast<char *>(framebuffer.data()), 3*height*width);
+    ofs.close();
 
     return 0;
 }