In previous experiments, we developed our 3D ray-tracing graphics engine to handle complex shapes such as diamonds. Those efforts begged an interesting question: what would it take to accurately simulate additional effects? In particular, I am interested in:

  • Transparent objects focussing light, in the way a magnifying glass can be used to create a bright point.
  • Diffuse reflections, in the way the underside of a ping-pong ball on a sheet of paper will be lit from beneath by light bouncing off the page.

Photons#

I can only thing of one way to do this: simulate individual photons in order to generate information about where they end up. In the art, this is known a photon map.

In this iteration of the ray-tracer, light sources do not illuminate objects directly. Instead, they emit a certain number of photons, with a certain amount of energy, in a certain pattern. We use the existing light simulation rules (absorption, reflection, refraction, or a combination of the three) to figure out where the energy the photon had ends up. We then need to process this set of points into an image. You might expand each point into a circle and blend neighbours together, or you might interpolate across a surface. Finally you can do your normal ray-tracing from the camera into the scene to determine what it all looks like.

With a photon map, light is used to 'paint' the surface of objects

The image above demonstrates the lensing effect (or caustic) of light being focussed by the glass sphere (centre). It also shows an interesting side-effect of my attempts to speed up computation by using Manhattan distance instead of Euclidian distance (saving millions of square roots) – objects are given a strange cross-hatching texture that gives the appearance that they are covered in canvas.

A Python module written in C++#

This version of the ray-tracer is built as a C++ module for Python. I have found this to be an excellent way of having your cake and eating it. Performance critical code can be ported to C++, whilst the user can continue to write Python to interact with it, even in a Jupyter notebook.

I used pybind for this. You end up with a bindings.cpp file which exposes various methods.

// bindings.cpp

#include <pybind11/pybind11.h>
#include <iostream>
#include <sstream>
#include <pybind11/iostream.h>
#include <pybind11/stl.h>

#include "core/Vector3D.h"
#include "core/Camera.h"
#include "core/Scene.h"
#include "objects/Object.h"
#include "objects/Sphere.h"
#include "core/Light.h"
#include "tools/Maths.h"
#include "tools/KDTree.h"
#include "tools/FLANN.h"
#include "objects/Material.h"


namespace py = pybind11;

PYBIND11_MODULE(rayTracerv2, m) {
    m.doc() = "Ray Tracer v2 module";

    py::add_ostream_redirect(m, "ostream_redirect");

    py::class_<Vector3D::Vector>(m, "Vector3D")
        .def(py::init<>())
        .def(py::init<float, float, float>())
        .def(py::init<bool>())
        .def("describe", &Vector3D::Vector::describe)
        .def_readwrite("x", &Vector3D::Vector::x)
        .def_readwrite("y", &Vector3D::Vector::y)
        .def_readwrite("z", &Vector3D::Vector::z);

    py::class_<Render3D::Camera>(m, "Camera3D")
        .def(py::init<>())
        .def(py::init<int, Vector3D::Vector, Vector3D::Vector, int, int, float>())
        .def("pre_calc", &Render3D::Camera::preCalc)
        .def("set_offset_grids", &Render3D::Camera::setOffsetGrids)
        .def("trace_rays", &Render3D::Camera::traceRays)
        .def_readwrite("id", &Render3D::Camera::id)
        .def_readwrite("position", &Render3D::Camera::position)
        .def_readwrite("direction", &Render3D::Camera::direction);

    py::class_<Render3D::Scene>(m, "Scene3D")
        .def(py::init<>())
        .def(py::init<int, std::vector<Render3D::Object*>>())
        .def_readonly("id", &Render3D::Scene::id)
        .def_readonly("n_objects", &Render3D::Scene::n_objects);

    py::class_<Render3D::Material>(m, "Material")
        .def(py::init<float, float, float, float>());

    py::class_<Render3D::Object>(m, "Object");

    py::class_<Render3D::Sphere, Render3D::Object>(m, "Sphere")
        .def(py::init<int, Vector3D::Vector, float, Render3D::Material*>());

    py::class_<Render3D::Cylinder, Render3D::Object>(m, "Cylinder")
        .def(py::init<int, Vector3D::Vector, Vector3D::Vector, float, Render3D::Material*>());

    py::class_<Render3D::Ray>(m, "Ray")
        .def(py::init<Vector3D::Vector, Vector3D::Vector, float>())
        .def("trace", &Render3D::Ray::trace);

    py::class_<Render3D::Photon>(m, "Photon")
        .def(py::init<>())
        .def_readonly("position", &Render3D::Photon::position)
        .def_readonly("energy", &Render3D::Photon::energy);

    py::class_<Render3D::PhotonMap>(m, "PhotonMap")
        .def(py::init<>())
        .def_readonly("n_photons", &Render3D::PhotonMap::n_photons)
        .def_readonly("photons", &Render3D::PhotonMap::photonTermini);
    
    py::class_<Render3D::PointLight>(m, "PointLight")
        .def(py::init<int, Vector3D::Vector, float, int>())
        .def("add_photons_to_map", &Render3D::PointLight::addPhotonsToMap);

    py::class_<KDTree::Tree>(m, "KDTree")
        .def(py::init<>())
        .def(py::init<const std::vector<Render3D::Photon>&>())
        .def("build_tree", &KDTree::Tree::buildTree)
        .def("range_search", &KDTree::Tree::rangeSearch)
        .def_readonly("root", &KDTree::Tree::root);

    py::class_<FLANN::PhotonMap>(m, "FLANNPhotonMap")
        .def(py::init<const std::vector<Render3D::Photon>&>());

    py::class_<KDTree::Node>(m, "KDNode")
        .def_readonly("point", &KDTree::Node::photon)
        .def_readonly("left", &KDTree::Node::left)
        .def_readonly("right", &KDTree::Node::right);

    m.def("rand_between", &Tools::randA, "random between -1 and 1");
    m.def("rand_vector", &Tools::flatRandVector, "random between -1 and 1");
}

Example Usage:# 

import numpy as np
import matplotlib.pyplot as plt

from rayTracerv2 import Vector3D, Ray, Camera3D, Scene3D, Material, Sphere, Cylinder, PhotonMap, PointLight, KDTree, KDNode

dusty_glass = Material(0.2, 0.2, 0.6, 1.5)
dusty_mirror = Material(0.2, 0.8, 0, 1)
dusty_metal = Material(0.8, 0.2, 0, 1)

sphere_1 = Sphere(0, Vector3D(0, 0.9, 0), 1.2, dusty_glass)
sphere_2 = Sphere(1, Vector3D(-2.5, 0.6, 0), 1.2, dusty_mirror)
sphere_3 = Sphere(2, Vector3D(2.5, 0.6, 0), 1.2, dusty_metal)

cyl_1 = Cylinder(3, Vector3D(0, -1, 0), Vector3D(0, -1.5, -0.1), 2, dusty_metal)
cyl_2 = Cylinder(4, Vector3D(-1.5, 1.3, 0), Vector3D(1.5, 0.2, 0), 1, dusty_glass)

scene = Scene3D(0, [cyl_1, sphere_1, sphere_2, sphere_3])

width, height = 400, 400
position = Vector3D(0, 0, 5)
direction = Vector3D(0, 0, 1)

camera = Camera3D(0, position, direction, width, height, 1.5)

camera.pre_calc()
camera.set_offset_grids()

PHOTON_RADIUS = 0.1

photon_map = PhotonMap()

point_light_a = PointLight(0, Vector3D(0, 14, 0), 1, 100_000)
point_light_a.add_photons_to_map(photon_map, scene)

photon_kdtree = KDTree(photon_map.photons)
pixels = camera.trace_rays(scene, photon_kdtree, PHOTON_RADIUS, 1000)

plt.figure(figsize=(10, 8))
plt.imshow(np.array(pixels).reshape((height, width)), cmap='grey')
plt.show()