This stuff is old.

Now that I've got this blog started, I thought it was time to retire my old homepage. It had a good innings (much of it dates from 1998) but it's time has come. The thing is, there were a few pages I had put up which I thought might just be worth dedicating a little space to here. Amongst these was a page for a simple raytracer I wrote back in 2001 and one for a basic radiosity renderer I wrote in 2006. So I spent a while trying to remember the bits and pieces I know about light/rendering and surfing various optics-related web pages and ended up with the below. I've tried to informally sketch the basics of global illumination and I hope it will be useful/interesting to my readers.

Models of light


Like many others, I can't resist indulging myself with some general (mostly irrelevant) waffle.

Classically (i.e. neglecting quantum effects) light is described by Maxwell's equations. Since I am a gauge theorist, I like to write these using differential forms: $$ \begin{align} dF &= 0 \quad\mbox{(this is just the Bianchi identity, there is no physics here)}\\ d*F &= 4\pi J \end{align} $$ Here F is the curvature of a U(1)-connection and is given by: $$ F = (E_xdx + E_ydy + E_zdz)\wedge cdt + (B_xdy\wedge dz + B_ydz\wedge dx + B_zdx\wedge dy) $$ where E and B are the usual electric and magnetic fields. J is the current 3-form and is given by: $$ J = -(j_xdy\wedge dz + j_ydz\wedge dx + j_zdx\wedge dy)\wedge cdt + \rho dx\wedge dy\wedge dz $$

In free space light is thus represented by the solutions to \( d*F = 0\) which is in fact the wave equation in disguise. As such, light is a wave. However in situations where the wavelength of the light is far smaller than the macroscopic dimensions and where the time scale we are interested in is much longer than the inverse frequency of the light (e.g. visible light in most everyday situations) a much simpler (and much older) model of light is applicable, namely geometric optics.

Geometric optics

In this model light is represented by rays which obey Fermat's principle of least time. It is possible to derive this model (via the eikonal equation) from Maxwell's equations and I hope to sketch this in a future post. For now I suggest we just accept geometric optics. It is unable to describe certain properties of light which are intrinsic to its wavelike nature (most notably diffraction and interference). This means that there are certain phenomena which this model cannot handle, e.g. thin film interference as seen in iridescence. Perhaps less familiar examples are Poisson's spot and Hamilton's celebrated conical refraction as exhibited by biaxial crystals such as aragonite or monoclinic double tungstate KGd(WO4)2.

Nevertheless, geometric optics captures most of our everyday experiences with visible light (indeed I think we our usually surprised by those phenomena which cannot be explained by geometric optics) and so we can expect to obtain reasonable images if we use it as the basis of our rendering.

Properties of light


The aptly named geometric optics talks only about the geometry of light rays. We must associate some physics to these rays in order to solve the rendering problem. The problem which the human eye/brain solve is to estimate the (time-averaged) amount of energy per second falling on each point of the retina as a function of wavelength (for those wavelengths in the visible spectrum). As such it is this physical quantity, energy per second, which we must associate to and calculate for our light rays.

Let us simplify things and ignore the lens and assume that the pupil is a point. Consider the light falling on a given point of the retina. Follow the ray of light which is illuminating this point backwards out into the scene which the eye is observing until we reach the source of the ray. This will be a point on some surface in the scene. The property of the surface which quantifies the rate at which energy is leaving it at a given point in a given direction is known as the radiance of the surface at this point in this direction. Note that we will also regard radiance as a property of this ray and so we will speak both of the radiance of a point of a surface in a given direction and also simply of the radiance of a ray of light in our scene. In any case, radiance is what we need to calculate. Evidently it has the dimensions of energy per unit time per unit area per steradian. (Incidentally, these units are equivalent to intensity per steradian but unfortunately the word intensity has been given similar but different meanings in different fields, including radiometry, and so it is best to avoid using it.)


So our retinas are able to estimate the energy spectrum of the light falling on them. Not bad. However just before we get too proud of ourselves it's worth noting that there is an important property of light that our eyes are almost completely insensitive to, namely the polarization of light (I say almost because of Haidinger's brush). For this reason almost all renderers ignore the polarization of light, thus greatly simplifying things. This does mean however that there are certain optical phenomena which we will never be able to render in our model. For example, we won't be able to render the way (appropriately oriented) polaroid sunglasses darken the sky or the light reflected from a flat surface like a calm lake much more than they darken the ambient light which has been scattered off diffuse surfaces. A more exotic example is birefringence as seen through calcite crystal. Nevertheless I think we'll get by!

Incidentally almost all animals, like humans, are essentially polarization-blind but there several striking exceptions. For example the eyes of cuttlefish are sensitive to polarization.

Cuttlefish would be much less satisfied with our decision to ignore polarization than we are.

Reflection and transmission of light


Fermat's principle tells us what happens when a ray of light strikes a surface: some of it will be reflected and some of it will be transmitted (i.e. refracted). The angle of reflectance is equal to the angle of incidence (Heron of Alexandria is credited with first noting this): $$ \theta_r = \theta_i $$ and angle of transmission is determined by the refractive index associated to the two media on either side of the surface, i.e. we have Snell's law: $$ \frac{\sin\theta_t}{\sin\theta_i} = \eta $$ What Fermat's principle does not tell us however, is how much of the light is reflected, how much is transmitted and in particular how these two quantities depend on the angle of incidence. It would be unkind to blame this on Fermat though; the ray model of light is completely geometric and so does not 'know' anything about the energy which its rays may be carrying. This means that in order to determine how much light is transmitted and reflected we have to briefly step outside our model, back to electromagnetism, and perform the calculation there. We can then bring the results home with us as some additional rules to use in our ray model of light.


The result we need is attributed to Fresnel. (Incidentally Poisson's spot was discovered as a result of Poisson's skepticism of Fresnel's findings in Fresnel's prize-winning 1818 essay on the wave nature of light.) Fresnel calculated the coefficients of reflection and transmission when the electromagnetic field meets a surface of discontinuity in the dielectric constant and obtained the Fresnel equations. These equations tell us the proportion of light reflected or transmitted as a function of the angle of incidence. There are two equations depending on whether the light is polarized perpendicular to or parallel to the plane of incidence (i.e. the plane in which the reflection/transmission occurs). For unpolarized light, we simply take the average of the coefficients for the two polarization states. More specifically if we let: $$ \begin{align} r_p(\theta_i) &= \frac{\eta\cos\theta_i - \cos\theta_t}{\eta\cos\theta_i + \cos\theta_t} \quad\mbox{(parallel polarization)}\\ r_s(\theta_i) &= \frac{\cos\theta_i - \eta\cos\theta_t}{\cos\theta_i + \eta\cos\theta_t} \quad\mbox{(perpendicular polarization)} \end{align} $$ then we coefficient of reflection we seek is given by: $$ F_r(\theta_i) = \frac{1}{2}\left(r_p(\theta_i)^2 + r_s(\theta_i)^2\right) $$ If Li is the radiance of the incoming ray, the radiance of the reflected ray is \( F_r(\theta_i)L_i\). Similarly, the coefficient of transmission is \( F_t = 1 - F_r\) (because of energy conservation) and the radiance of the transmitted ray is \( F_t(\theta_i)L_i\).

Note that the above equations are appropriate for dielectric materials (i.e. materials which do not absorb light, e.g. perfect glass). For conductors (i.e. materials which do absorb light, e.g. metal) the index of refraction acquires an imaginary component. The above equations can still be used (and Snell's law still holds, even with the non-real index of refraction) but some extra care is needed. In any case, the point is that there is no problem dealing with such materials.


'Great,' I hear you say, 'thanks to Fresnel and Fermat we can calculate all we need to know when a ray of light strikes a surface.' It is true that if we follow the above we can calculate what happens when light meets perfectly smooth surfaces like mirrors and blocks of glass but it would be nice to be able to render more than just mirrors and blocks of glass. The sort of reflection and transmission which occurs when light meets mirrors and glass is known as perfect specular reflection and transmission. We want to be able to handle more general types of light reflection (from now on we'll just talk about reflection, it is not really any harder to handle transmission).

There are lots of reasons why most surfaces do not exhibit perfect specular reflection. For example the surface may be rough and so, although on average the normals will follow the macroscopic shape of the surface, they will be randomly distributed (according to some probability distribution). When light strikes such a surface it may in principle be reflected in any direction. As a result the coefficient of reflection Fr no longer depends only on the angle of incidence, it also depends on the angle of reflection. Furthermore the reflected ray may not lie in the plane of incidence and so the amount of light reflected can also depend on the angle \( \phi\) between the plane of incidence and the plane of reflection (i.e. the plane containing the reflected ray and the surface normal). (Experts will note that I'm excluding anisotropic surfaces.) In this setting the coefficient of reflection is usually denoted in lowercase as fr and is known as the bidirectional reflectance distribution function or BRDF for short. Using the BRDF, if we know that a ray of light with radiance Li hits a surface with angle of incidence \( \theta_i\) and leaves with angle of reflection \( \theta_r\) and that there is an angle \( \phi\) between the planes of incidence and reflection, then the radiance of the reflected ray is: $$ L_r = f_r(\theta_i, \theta_r, \phi)L_i $$

We can express what we know about perfect specular reflection using a BRDF in terms of the Dirac delta distribution: $$ f_r(\theta_i, \theta_r, \phi) = \frac{\delta(\theta_i - \theta_r)\delta(\phi)F_r(\theta_i)}{\cos\theta_i} $$ Evidently we need BRDFs for surfaces which are not perfect specular reflectors. Fortunately there is no shortage of these. Firstly, at the opposite extreme to the perfect specular BRDF, there is the perfect diffuse BRDF: $$ f_r(\theta_i, \theta_r, \phi) = \mbox{const.}\cos\theta_r $$ We should probably also mention the Blinn-Phong 'BRDF' (which is not actually a BRDF since it is not 'bidirectional') which underlies the standard lighting models used on modern graphics hardware): $$ f_r(\theta_i, \theta_r, \phi) = \mbox{const.}\frac{\cos(\alpha)^n}{\cos\theta_i} $$ (where \( \alpha\) is the angle that the ray lying halfway between the incident and reflected rays makes with the normal.) After these come the more realistic BRDFs which can still be expressed using analytical formulae. These can be phenomenological e.g. the Schlick BRDF or they can be the result of modelling what is happening for surfaces that are not perfect specular reflectors. For example, we can model a rough a surface by imagining that the surface is really made up of a many small perfect perfect specular reflectors which are randomly oriented according to some chosen probability distribution. Such models are often called microfacet surface models and include the Torrance-Sparrow BRDF (introduced in the field of heat transfer in 1967) and the Cook-Torrance BRDF amongst many others. Another example of a well known BRDF is the Oren-Nayar BRDF in which a Gaussian distribution of perfectly diffuse microfacet reflectors is assumed. A relatively simple formula gives a good approximation to the resulting BRDF so I might as well give it. It depends only on the variance \( \sigma^2\) of the Gaussian distribution (note that \( \sigma\) is measured in radians). $$ f_r(\theta_i, \theta_r, \phi) = const.\left(A + B \max(0, \cos\phi\sin(\max(\theta_i, \theta_r))\tan(\min(\theta_i, \theta_r))\right) $$ where: $$ \begin{align} A &= 1 - \frac{\sigma^2}{2(\sigma^2 + 0.33)}\\ B &= \frac{0.45\sigma^2}{\sigma^2 + 0.09} \end{align} $$ In addition there are of course empirical BRDFs. These days there are free online databases of BRDFs which are the result of sampling the BRDF of real world substances like plastics, paints, metals, fabrics etc. For example there is Cornell's database, the Columbia-Utrecht Reflectance and Texture Database, the Bonn BTF Database and Wojciech Matusik's data.

We should also say that the fact that a surface is not completely smooth is far from the only reason it will fail to exhibit perfect specular reflection. In general, when light strikes a surface it will enter the surface and exit at a different point. It is also possible to model this subsurface scattering phenomenon but it is significantly more costly to compute. Translucent substances like skin and marble are typical examples of materials which exhibit significant subsurface scattering.

Raytracing and the rendering equation

Imagine holding up a rectangular pane of glass in front of your eyes and observing a scene through it. In fact let's simplify and assume our eyes are a single point. Then consider all the rays of light in the scene arriving at this point through the glass. We want to calculate the radiance of each of these rays. Once we have computed this, we will ask our computer to light up the pixels on its screen according to these radiance values. The light emitted by the computer screen will thus be the same as the light that would pass through the pane of glass if we were really observing the scene with our computer screen playing the role of the pane of glass.

In raytracing we follow each ray backwards from the point representing our eyes through the screen till we hit a point on a surface in the scene. We must calculate the outgoing radiance Lo of this point on the surface. This is a sum of two terms: $$ L_o = L_e + L_r $$ Here Le is the radiance of any light which the surface may be emitting (i.e. if it is a light source) and Lr is the radiance of any light which the surface may be reflecting at this point in the outgoing direction. As noted above, the reflected light can be the result of incident light from any direction and so to calculate it we must add up, i.e. integrate, the reflected light over the hemi-sphere of all incident directions. We thus have: $$ L_r(\theta_r) = \int_0^{2\pi}d\phi\int_0^{\pi/2}\sin\theta_i d\theta_i f_r(\theta_i, \theta_r, \phi)L_i(\theta_i, \phi) \cos\theta_i $$ This integral is usually estimated using Monte Carlo methods (a notable exception occurs in radiosity rendering). In order to do this we need to sample (i.e. calculate) \( L_i(\theta_i, \phi)\). Since the incident light on one surface is just the outgoing light for another, we calculate Li just as we calculate Lo above. Thus for a given direction specified by \( \theta_i, \phi\) we follow the ray from the point of our surface off in this direction till we hit another point of a surface in the scene and calculate the outgoing light in this direction using the above. Obviously this is a recursive definition and in order to avoid perpetually following reflected light rays around the scene, we must specify a rule for when we stop tracing reflected light. In general this is a difficult problem because parts of the scene will look fine when rendered with only a small number of bounces whereas other regions will need a great many. Perhaps the simplest rule (which I used in my raytracer below) is to set a maximum number of bounces for all rays and stop when this is reached.

The above system for calculating radiance can be neatly summarised using the rendering equation. The rendering equation was simultaneously introduced into computer graphics by Kajiya and Immel et al. at SIGGRAPH 1986. The rendering equation is essentially the combination of the equations for Lr and Lo above, except that instead of integrating over an abstract hemi-sphere of directions, the integral is taken over the points of the surfaces in the scene. In other words, we can parameterise the rays incident on a point of a surface either by their direction as above or by the point of the surface from which they originate. (Incidentally, the above equations for Lr and Lo are sometimes referred to as the rendering equation. In my opinion it is worth reserving the term rendering equation for the equation in which the integral is taken over the surfaces in the scene.) The only thing we need to do then is to express the area element dS of a surface in the scene in terms of the solid angle form \( \sin\theta_i d\theta_i\wedge d\phi\) used above. This is: $$ \sin\theta_i d\theta_i\wedge d\phi = \frac{\cos\theta'_o}{r^2}dS $$ where \( \theta'_o\) is the angle that the ray illuminating our point makes with the normal of the surface from which it is an outgoing ray. For convenience (and following Kajiya) we introduce a geometry term: $$ G = \frac{\cos\theta'_o \cos\theta_i}{r^2} $$ and a visibility term: $$ V(x, x') = \begin{cases} 0 & \mbox{if the points $x$ and $x'$ are ocluded by the geometry}\\ 1 & \mbox{otherwise} \end{cases} $$ Writing our previous equations in this notation gives the rendering equation: $$ L_o = L_e + \int_S f_rVGL_o $$ (assuming the distribution of light in our scene has reached steady state.) One of the nice things about formulating things this way is that if we define the integral operator T by: $$ TL = \int_S f_rVGL $$ then we can formally manipulate the rendering equation as follows: $$ L_o = L_e + TL_o\\ (I - T)L_o = L_e\\ L_o = (I - T)^{-1}L_e = (I + T + T^2 + \cdots) L_e $$ The above expansion in powers of T has the physical interpretation that the outgoing light at a point is the sum of the emitted light which as reached that point after 0 bounces (i.e. the light that point itself emits) plus the emitted light which has reached that point after one bounce, two bounces etc. In other words, the power of T counts the number of times that the emitted light has been reflected.

Having finally reached the rendering equation, it's seems appropriate to present some of my own efforts at solving it. Unfortunately, writing a full-blown Monte Carlo raytracer remains on my TODO list but as I mentioned I wrote two simple renderers some years ago which yield approximate solutions.

The simple raytracing approximation

I have already mentioned that I wrote a simple raytracer in 2001. There is of course no shortage of simple raytracer source code in the world but I still thought it would do no harm to make the source code to mine available. You can find it here. It is written in C++, is self contained and should work with almost any compiler/platform combination; it was written using GNU C++ under Linux (unfortunately I don't know the compiler version but it was probably the newest in 2001 when I wrote this).

This is a very simple type of raytracer. The BRDF for each surface is a sum of perfect diffuse component term as well as a perfect specular Phong component. Specular light is reflected correctly around the scene but only direct diffuse light is calculated (i.e. diffuse interreflections are ignored but specular interreflections and diffuse to specular reflections are handled correctly). A constant 'ambient' light term is also introduced to take some account of the missing interreflected diffuse light. Although the lighting does not look physically realistic I was still quite pleased with the results at the time I wrote this.

As is evident, each image is a view of the same scene, looking from a different location. The characteristic hard shadows and perfect specular reflections of this type of raytracing are clearly visible.

The radiosity approximation

In radiosity rendering all surfaces are assumed to be perfectly diffuse. The surfaces in the scene are subdivided into a collection of small surface elements which are each assumed to have constant radiance. In this case the rendering equation becomes finite dimensional and so in principle it can be solved exactly by calculating the matrix of T (above). Furthermore because we have a perfectly diffuse BRDF, this matrix depends only on the geometry of the surface elements, and so in particular it does not depend on the location of the viewer. For this reason we can quickly rerender the scene looking from a different location or in a different direction (this is why lightmaps work).

I have mentioned that to solve the rendering equation in this setting, we merely need to compute and invert a matrix. The rows and columns of this marix are indexed by the surface elements in the scene. An entry in this matrix (called a form factor) is thus determined by a pair of surface elements and corresponds to the geometric term which quantifies the extent to which they exchange diffuse light with each other. Since the vast majority of pairs of surface elements have little effect on each other, it is an inefficient approximation to calculate this full matrix of form factors. This was recognised soon after radiosity was introduced and progressive radiosity (using the hemi-cube method was introduced). Using these techniques we spend most of our time calculating the most important form factors, we do this efficiently using rasteriser triangle filling techniques and we avoid inverting the matrix by using the Gauss-Seidel method.

A good place to read about all of the above, i.e. radiosity, the hemi-cube method and the progressive radiosity technique is the original papers:

As I mentioned, I wrote a simple radiosity renderer in 2006. Again, I thought it would do no harm to make the source code available here (note that the simple 3DS loader included is not my work, it was one I found on the web and appears to have written by a chap called Jerry Chen). I wrote it in C++ using Visual Studio C++ 8.0 2005 Express Edition. It needs OpenGL (note that OpenGL is not used to do any lighting, just to plot triangles).

Some sample images I generated are below. Note that because diffuse light is handled correctly, the shadows look far more realistic than for my raytracer. Since progressive radiosity produces images which converge to the final solution, it is possible to watch the algorithm as it converges. The four images below show this process with the fourth image being essentially converged. (Some aliasing problems are visible.)

The next image below shows how the fourth image above is actually represented as a collection of constant colour triangles (recall the radiosity assumption is that surface elements each have constant radiance). The smooth images are generated by using linear interpolation.

Taking advantage of the fact that the radiosity solution can be easily rerendered from a different location in the scene, I also rendered the solution displayed in the fourth image above from two different locations.

The below two images are the same scene with brighter lights. I like the shadows in these images.

This last image is the same scene with an extra yellow light shining.