A JavaScript RayTracer

If you don’t want to hear about all this and go straight to the raytracer, you can find it here: live demo.

Models of Light

Light is a very complex system to model accurately, which is why it is so difficult to generate photorealistic images on a computer. As is always the case, more complex and realistic simulations require more computing power and therefore reduce speed at which the simulation will run. Computer scientist have been working hard for a long time to improve on the performance of algorithms and find good payoffs between model accuracy and computational complexity.

Before we dive into the technicalities of computer graphics in general and raytracing in particular, let us take a moment and try to understand some things about light.

Maxwellian description

The Maxwell’s equations are a set of partial differential equations that, together with the Lorentz force law, form the foundation of classical electrodynamics:

(1)   \begin{align*} \nabla \vec{E} &= \frac{\rho}{\epsilon_0}, \quad \nabla \times \vec{E} = -\frac{\partial \vec{B}}{\partial t} \nonumber \\ \nabla \vec{B} &= 0, \quad \nabla \times \vec{B} = \mu_0 \left( \vec{j} + \epsilon_0 \frac{\partial \vec{E}}{\partial t} \right) \end{align*}

In regions where no charges are present (\rho = 0) and no currents exist (\vec{j} = 0), these can be brought into the form of two wave equations for the electric and magnetic field:

(2)   \begin{equation*} \frac{1}{c^2} \frac{\partial^2 \vec{E}}{\partial t^2} - \nabla^2 \vec{E} = 0, \quad \frac{1}{c^2} \frac{\partial^2 \vec{B}}{\partial t^2} - \nabla^2 \vec{B} = 0 \end{equation*}

This is a second order partial differential equation that describes the propagation of electromagnetic waves through the vacuum with c the speed of propagation, explicitly given by (for the vacuum case)

(3)   \begin{equation*} c = \frac{1}{\sqrt{\epsilon_0 \mu_0}}, \end{equation*}

where \epsilon_0 is the vacuum permittivity and \mu_0 the vacuum permeability. For the simplest model of real medium the wave propagates in, \epsilon_0 and \mu_0 are multiplied by material specific constants \epsilon_r and \mu_r, which determine the optical density of the medium (see below).

Maxwell’s classical theory therefore describes light as waves which propagate through media subject to the differential equation above. However this (microscopic) level of description is hardly ever used in computer graphics. Here one is usually in the situation where the wavelength of light is small compared to the size of the structures with which the light interacts. Consequently, for most imaging applications it is sufficient to describe the propagation of light on the level of geometrical optics.

Geometrical optics

In geometrical optics light is not described as waves but rather as rays, which propagate (in homogenous media) in straight lines. The ray here is an idealisation used to model how wavefronts of light would propagate through space. The model relies on Fermat’s principle (a consequence of Huygens’ principle), which can be stated in its modern formulation as

The optical length of the path followed by light between two fixed points, A and B, is an extremum. The optical length is defined as the physical length multiplied by the refractive index of the material.

As such light rays are defined to propagate in a rectilinear path as they travel in a homogeneous medium. Rays bend (and may split in two) at the interface between two dissimilar media, may curve in a medium where the refractive index changes, and may be absorbed and reflected. Thus, compared to Maxwell’s equations the level of complexity for the description of light has greatly reduced: instead of solving a complicated boundary value problem, one reverts to the reconstruction of the most important pathways. For this in principle a calculator and a ruler suffice.

This simplification however comes at a price: Essentially the wave-like nature of light has been discarded and as such certain effects, e.g. diffraction and interference, are not accounted for anymore. The following beautiful example of iridescence cannot be understood on the level of geometrical optics thereby.

iridiscence

Perception of light

The human eye is essentially a spectrum analyser as it estimates the frequency-resolved, time-averaged amount of energy per unit time that falls onto a given point of the retina. Under the simplifying assumption that we model the eye as a zero-dimensional point, we can associate with each ray that hits the eye a direction. Following the ray backwards out into the scene we reach a point on a source surface that the ray was emitted from. The physical \textit{property of the surface} that determines the quantity of radiation that is emitted from a given point and falls into a solid angle in a specified direction is termed radiance. The SI unit of radiance is watts per steradian per square metre (\mathrm{W~sr^{-1}~m^{-2}}). As such, the task of any computer graphics system is to calculate this physical quantity for all visible points of the surfaces in the scene.

Before we move on it should be mentioned that there is another important property of light which, however, our human eye is blind to: polarisation. Electromagnetic waves propagating in three-dimensional space can oscillate in more than one direction: after fixing the direction of propagation the electro-magnetic fields can oscillate in any one direction contained in the two-dimensional plane perpendicular to this direction (note that sound waves for instance have only one sense of polarisation). By convention, the polarisation of light refers to the polarisation of the \textit{electric} field. In general the solution of the wave equation for a wave propagating along the z-direction can be written as

(4)   \begin{align*} \mathbf{E}(\mathbf{r},t) &= |\mathbf{E}| \text{Re} \left\{ \vec{\Psi} \exp(i (kz - \omega t)) \right\} \nonumber \\ \mathbf{B}(\mathbf{r},t) &= \mathbf{\hat{z}} \times \mathbf{E}(\mathbf{r},t) \end{align*}

with \omega = c k the angular frequency of the wave. \vec{\Psi} is the Jones vector given by

(5)   \begin{equation*} \vec{\Psi} = \begin{pmatrix}\cos(\theta) \exp(i \alpha_x)\\\sin(\theta) \exp(i \alpha_y)\end{pmatrix} \end{equation*}

For \alpha_x = \alpha_y the fields oscillate along a single direction (determined by \theta) commonly referred to as linear polarisation. For \alpha_x = \alpha_y + \frac{\pi}{2} the field vector rotates at the optical frequency \omega in a circle which is called circular polarisation.

Although the human eye is not sensitive to the polarisation of light, there are animals which have evolved eyes that can determine the sense of polarisation. A famous example is the honeybee (but also other invertebrates): bees communicate the location of a source of food by means of a dance that points in the direction of the goal. Considerable distances to the food suggested that bees do not fix their bearings on the local terrain. Indeed it was found that blocking a bees view of the sky prevents their proper orientation. On the other hand if the bees are given as little as one percent view of the sky, they are able to communicate the location of the food accurately. The resolution to this puzzling observation is the fact that the earth’s atmosphere acts like a spatially dependent polarisation filter (which creates a series of concentric ring arcs in the celestial hemisphere that are parallel to the orientation of the skylight’s \mathbf{E} vector and perpendicular to the great circles that pass through the sun) and bees can detect the pattern of ultraviolet , polarised light in the sky.

BEES_1

Typically the output of a computer graphics system is intended for consumption by human eyes. Therefore in the following we will make the simplifying assumption that all light sources emit unpolarised light and all scattering of light from any surface preserves this property.

Reflection and transmission of light

When light hits a surface it is partly reflected and partly transmitted. The simplified propagation model of geometrical optics accounts for this observation by two rays which leave the intersection point: (a) a reflected ray, which is emitted from the point of reflectance at an angle \theta_r equal to the angle of incidence \theta_i

(4)   \begin{equation*} \theta_r = \theta_i \end{equation*}

(b) a transmitted ray which leaves the intersection point at an angle \theta_t given by Snell’s law

(5)   \begin{equation*} \frac{\sin(\theta_t)}{\sin(\theta_i)} = \eta \end{equation*}

where \eta = \frac{n_i}{n_t} and n_i is the optical density of the medium the incident ray propagates in and n_t the optical density for the transmitted ray. Both of these follow from the application of Fermat’s principle. The information not provided by this principle is, however, the amount of light which is reflected.

Snells_law

Fresnel

At the beginning of the 19th century French engineer and physicist Augustin-Jean Fresnel studied the behaviour of light both theoretically and experimentally. His findings are summarised in the \textit{Fresnel equations} which describe the fraction of the incident power that is reflected from (R) and transmitted through (T) an interface between two media of different optical density as a function of the angle of incidence. The calculations of R and T depend on polarisation of the incident ray. Two cases are analyzed:

  • The incident light is polarized with its electric field perpendicular to the plane containing the incident, reflected, and refracted rays. In this case the light is said to be s-polarized.
  • The incident light is polarized with its electric field parallel to the plane described above. Such light is described as p-polarized.

Fresnel found that the reflectance is given by

(6)   \begin{align*} R_s(\theta_i) &= \left| \frac{n_1 \cos \theta_i - n_2 \sqrt{1- \left( \frac{n_1}{n_2} \sin \theta_i \right)^2}}{n_1 \cos \theta_i + n_2 \sqrt{1- \left( \frac{n_1}{n_2} \sin \theta_i \right)^2}} \right|^2 \nonumber \\ R_p(\theta_i) &= \left| \frac{n_1 \sqrt{1- \left( \frac{n_1}{n_2} \sin \theta_i \right)^2} - n_2 \cos \theta_i}{n_1 \sqrt{1- \left( \frac{n_1}{n_2} \sin \theta_i \right)^2} + n_2 \cos \theta_i} \right|^2 \end{align*}

and the transmission coefficients obey (due to energy conservation) T_s = 1 - R_s and T_p = 1 - R_p. For unpolarised light one simply takes the average R = (R_s + R_p) / 2. For an incoming ray of radiance L_i that intersects the interface at an angle \theta_i the radiance of the reflected ray (which leaves the point of intersection and the angle \theta_t = \theta_i, see above) is given by R L_i.

Limitations of the Fresnel equations

The limitations of the Fresnel treatment of reflection and transmission are best appreciated when formulated in a much more general language. Enter \textit{bidirectional reflectance distribution functions} (BRDF’s). Imagine light from an infinitely far away source hitting a surface from direction \mathbf{\omega}_i (with \theta_i the angle between \mathbf{\omega}_i and the surface normal \mathbf{n}). The BRDF f(\omega_i, \omega_r) describes the ratio of reflected radiance exiting along a direction \omega_r to the irradiance incident on the surface from direction \omega_i

(7)   \begin{equation*} f(\omega_i, \omega_r) = \frac{d L_r(\omega_r)}{L_i(\omega_i) \cos(\theta_i) d\omega_i} \end{equation*}

with L_r the radiance reflected in direction \omega_r, L_i the incoming radiance from direction \omega_i. For isotropic surfaces (surfaces where the BRDF only depends on the inclination at which the incoming ray hits the surface and not the polar direction) this function can be parametrized by the angles (\theta_i, \theta_r the angle between \mathbf{\omega}_r and the surface normal \mathbf{n} and a polar angle \phi which describes the polar deviation between the outgoing direction \mathbf{\omega}_r and the plane of incidence, spanned by \mathbf{\omega}_i and n.

The BRDF that corresponds to the Fresnel theory is therefore given by

(8)   \begin{equation*} f(\theta_i, \theta_r, \phi) = \frac{\delta(\theta_i - \theta_r) \delta(\phi) R(\theta_i)}{\cos \theta_i}. \end{equation*}

For perfectly smooth, very shiny surface the above equation provides a sufficient idealisation, but most surfaces that surround us, cannot be described by such a simple BRDF. For rough surfaces one can think of a random distribution of normal vectors (with an average normal direction which coincides with the “macroscopic” normal direction). Therefore when the Fresnel BRDF is averaged over tiny portions of the surface, the effective BRDF will not be sharply peaked about the mirror-reflectance direction but washed out.

Another important effect is subsurface scattering (SSS). Here light penetrates the surface of a (partly) translucent object and is scattered within the material through interaction and then exits the surface at a different point. Subsurface scattering is important for the realistic rendering of materials such as marble, skin, leaves, wax and milk for instance.

Perfectly diffuse surfaces are described by the Lambertian reflectance BRDF which is simply given by f(\theta_i, \theta_r, \phi) \propto \cos \theta_r. The resulting distribution of outgoing radiance is therefore perfectly the same in all directions: once an incoming ray is reflected by a Lambertian reflector all information about its origin is completely lost.

Another example, especially popular in computer graphics is the Blinn-Phong shading model. Here the reflection of light from a surface is modelled as the combination of different effects. There is an ambient component which is independent from the direction of incidence, reflectance and also from the radiance of the incoming beam. The ambient component is simply given by L_{\rm ambient} = L_{\rm surrounding} R_{a}, where L_{\rm surrounding} is a constant radiance emitted by the surrounding (with no associated direction of origin) an R_a a reflection constant averaged over all surfaces.

A second contribution comes from diffuse reflection. Here light coming from point sources is reflected by what is assumed to be a perfect Lambertian scatter. Therefore the resulting reflected radiance does not depend the angle of observation \theta_r. It does however depend on the angle of incidence \theta_i and therefore on the position of the light relative to the reflecting surface as the illumination of the surface changes as a function of the angle of incidence. The reflected radiance is given by

(9)   \begin{equation*} L_{\rm diffuse} = L_{\rm in} R_{\rm diffuse} \mathbf{L} \cdot \mathbf{N} \end{equation*}

where L_{\in} is the radiance of the point source, R_{\rm diffuse} a surface material specific reflection constant, \mathbf{N} the surface normal and \mathbf{L} a unit vector point in the (negative) direction of the incoming ray (the illumination of the surface and therefore the diffuse component peaks when light hits perpendicularly).

The last contribution comes from perfect specular reflection in the sense of Fresnel. Here the shininess of the surface is accounted for. Here radiance of the reflected light depends on the direction of incidence, a surface material depend specular reflectance, and the direction of observation.

(10)   \begin{equation*} L_{\rm specular} = L_{\rm in} R_{\rm specular} (\mathbf{R} \cdot \mathbf{D})^n \end{equation*}

with \mathbf{R} the direction of ideal specular reflection and \mathbf{D} = \omega_r the direction of observation. The exponent n models the shininess of the surface.

Raytracing

Raytracing is an elegant technique which allows to determine the radiance of light inciding on the eye of the observer from a particular direction by following the ray backwards out into the scene and examining the surface it was emitted from. The goal is to calculate the outgoing radiance L_{\rm out} from that point in the particular direction. L_{\rm out} is the sum of two terms

(11)   \begin{equation*} L_{\rm out} = L_{\rm e} + L_{\rm r} \end{equation*}

where L_{\rm e} is the radiance of light that the surface itself is emitting (think about a lamp) and L_{\rm r} the radiance of the light incoming from other points in the scene and reflected by the surface. In order to calculate L_{\rm r} we integrate the BRDF over all incoming directions

(12)   \begin{equation*} L_{\rm r}(\omega_r) = \int_{\mathcal{S}^2} d\Omega_i f(\omega_i, \omega_r) L_i(\omega_i). \end{equation*}

Here of course L_i(\omega_i) being the incoming light in direction \omega_i is just the outgoing light L_{\rm out} of another surface. We may determine this quantity by following back from the point of incidence in the direction \omega_i and thus iterating the above equation. This is a recursive definition and may calculate an approximate solution by limiting the maximum recursion depth.

In the rendering equation the integral over all incoming directions is replace by an integral over all surfaces in the scene as those are the only direction light might be coming from.

Implementation

I have implemented the simplest version of these ideas in a ray tracer written in JavaScript. Ray tracing is quite a computationally expensive algorithm and therefore implementing it in a scripting language is a spectacularly bad idea. Nonetheless, the attractive aspects about it where that it can be directly delivered into the browser.

The BRDF used in the implementation is based on the Phong reflection model. It therefore comprises

  1. ambient contribution: modelling the surrounding isotropic illumination of the scene
  2. local diffuse reflection: modelling diffuse reflection of light emitted from point lights distributed in the scene
  3. local specular reflection: modelling light specular reflection of light originating from point lights
  4. reflected specular component: only the specular (mirror-like) reflection of light originating form other surfaces is taken into account. This contribution is calculated by recursively backtracking the ray along the specular direction.

Note that the first three components only reflect properties of the local surface and the point light sources. Only due to the last component properties of other surfaces (the colour for instance) contribute to the emitted light.

Rendering results

A live demo can be found here.

render