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.

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.

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.)

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.

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.

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 F_{r} 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 f_{r} 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 L_{i} 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.

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 L_{o} of this point on the surface. This is a sum of two terms:
$$
L_o = L_e + L_r
$$
Here L_{e} is the radiance of any light which the surface may be emitting (i.e. if it is a light source) and L_{r} 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 L_{i} just as we calculate L_{o} 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 L_{r} and L_{o} 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 L_{r} and L_{o} 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.

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.

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:

- Modeling the Interaction of Light Between Diffuse Surfaces
- The hemi-cube: a radiosity solution for complex environments
- A Progressive Refinement Approach to Fast Radiosity Image Generation

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.