Rays: Casting and Paraxial Optics

1. Paraxial Geometrical Optics and the System Matrix

Mathematical Models

"Now consider a spherical chicken..."

-- hoary punch line

This material is from the first four (out of 25) Chapters in [2], a typical sophomore-level textbook.

Light is an electromagnetic wave with many fascinating properties. Like also being an elementary particle. So much for Chapters 1 and 2 of [2].

Refraction is the bending of light by an optical medium, reflection is the bouncing-off of light from a surface. Geometrical optics ([2], Ch.3) is a special case of physical optics, which tosses overboard all the wave (and quantum in general) attributes of light. Christiaan Huygens (Dutch, b. 1629) had an ether-based theory of light that was a sort of geometric optical approximation in that it failed for diffraction but did derive the law of refraction. Pierre Fermat (French, b. 1601) used a minimum-principle argument (minimizing light's travel time) to prove the same law. His principle is also not quite right, since there turn out to be time-maximizing light paths and cases where all paths are equivalent in time (between the foci of an ellipsoidal mirror, for instance).

The law of refraction is named, for no very good reason, Snell's law. For one thing, his name wasn't Snell: it was Willebrord Snellius (Dutch, b. 1580). For another, the law as been known a long time, being investigated by Ptolemy (Greek, b. ca. 90 AD), and Witelo in the middle ages, but they didn't have the formalism to write it down and used tables instead. (In fact, the first human ancestor whose spear went over the fish knew about it). One Ibn Sahl (Persia, b. ca. 940 AD) discovered it (his manuscript is the frontispiece of this project page.) What is it, anyway?

The law of refraction (Snell's Law): When a ray of light is refracted at an interface between two uniform media, the transmitted ray remains in the plane of incidence and the sine of the angle of refraction is directly proportional to the sine of the angle of incidence.

There's also a law of reflection: When a ray of light is reflected at an interface dividing two uniform media, the reflected ray remains in the plane of incidence, and the angle of reflection equals the angle of incidence. The plane of incidence is defined by the incident ray and the surface normal vector at the point of incidence.

Snell's Law is usually written:
n1 sin θ1 = n2 sin θ2,
with n1, n1 the refractive indices of the two media, where the speed of light vin a medium of refractive index n is v = c/n if c is the speed of light in a vacuum. Slower lightspeed in second medium, more refraction. Here's a picture of reflection and refraction.


Now using these two laws we can analyze systems of lenses and mirrors like the all-lens system below: snell

After all this, we'd like to use Snell's Law in our geometrical optics approximation to reality. But we have a problem. Those sines are nonlinear functions. We won't be able to find solutions by the simple-minded linear methods we know and love. Well, then! Let's toss some more reality into the bin. Let's just consider only paraxial rays, those that never make large angles with the optical axis. In fact, those for which we can approximate,
with angles measured in radians, not degrees (!!)
sinθ = tan θ = θ
That sure cleans things up: now the paraxial Snell's Law is
n1 θ1 = n2 θ2 .

The paraxial assumption is an example of the "small angle approximation" so dear to physics and engineering. It is a linear approximation, which describes a given function (in some locality) as a linear function. Here, for example sin(x) &sin x near 0. In the future you'll often expand functions in infinite (Taylor) series and ignore all but the first term to get your linearization. It happens that sin(x) = x - x3/3! + x5/5! - ..., near 0, and for small x we're only considering the first term.

Now notice in the figure above that different rays are going to travel different distances through the lenses, which are fatter and thinner. So just what happens depends on just where and at what angle the ray enters the lens. Well, that's inconvenient. But wait! We know what to do! Just jettison that bit of reality. The thin lens assumption is: "consider an infinitely thin lens" with all the refractive power and none of the annoying size of an actual object. A diagram might look like this, with X for convex and V for concave lenses: thinlens

Using pretty simple algebra (but pretty complicated diagrams) [2] derives elegant formulae for useful, practical numbers that help in understanding and creating lenses. They're derived from the basic geometry (say radii of curvature, distance) and basic physics (snell's law, indices of refraction). The result are usable parameters that describe spherical mirrors, refraction at spherical surfaces, thin lenses etc.

An important thin lens concept is its focal length f, the (signed) distance to the image it forms of an object at infinity. Convex lenses have positive focal lengths, Concave negative). The lensmaker's equation predicts the focal length of a lens in terms of its refractive index, that of the medium it is in, and its radii of curvature. In air (of refractive index n = 1), one version of the thin-lens equation states that an object at distance s is imaged at point s' by a lens of focal length f if 1/s + 1/s' = 1/f. The power P of a lens measures how strongly it bends light, and is defined as 1/f. In an x-y plane, We measure angles (in radians) that rays deviate from the optical system (x) axis (rays parallel to the x-axis are at 0 radians). Rays from the left rising in the y direction are positive. A lens bends an incoming ray of angle θ, according to our approximations, to give it a new angle:
θ' = θ - yP
if the ray impinges on the lens at height y. Thus
θ' = θ -(1/f)y

That's our quickie summary of [2], Ch. 3. Remember angles are always in radians, not degrees. So now we can use our paraxial geometrical optical spherical chicken and do some (approximate) reasoning. Snideness aside, creating solvable linear approximations is a crucial skill you'll be developing in the next three years.

The Method

Often linear systems are described by the combination (summation, multiplication, convolution,..., or other linear operators) of transformations produced by the system's parts. We shall describe a ray moving left to right by its instantaneous height and angle of flight: [y,θ]T. We'll see that paraxial optical elements can be described by 2x2 matrix transformations (ray transfer functions) representing what happens to a ray at an optical element, and the product of the transfer matrices describes the optical system as a system transfer function. This material is covered in [2] Ch. 4, for example.

Diagram the system as above, a 2-D plot of Y versus X, with the X axis being the optical axis. An object being imaged is considered to be in an input plane at the left of the diagram (below, it's at x = 0), and the output plane is wherever we want to compute the final height of the ray. Things we might want to know about the system include (there are many more) its front and back focal length. The latter is just the effective focal length of the entire system, and the former is similarly the focus point of angle-0 rays coming in from the right, headed left. The axial image point of an object point nearer than infinity on the optic axis is where on the optic axis all its rays are focussed (yes, they all go through one point). The image plane is the plane (line in a 2-D diagram) through the axial image point perpendicular to the optic axis.

The linear magnification of the system is yf/y0, the relative height of a ray at the input and output planes.

Consider a ray [y0, θ0 ]T moving (translating) through a homogeneous medium. Angles are measured in radians, remember. If it moves axially by L at an angle θ0, Its new description is
[y0 + L tan θ0, θ0]T. The paraxial assumption is that tan(x) = sin(x) = x, so we can describe the resulting linear transformation of the ray's description by

|y1| = |1 L|  |y0|
|θ1| = |0 1|  |θ0|
This simple matrix is how we describe the change in the ray as if moves through some uniform medium for an axial distance.

Remember this equation from the end of the last section?
θ' = θ -(1/f)y.
So for a thin lens, guess what?

|y1| = |1    0|  |y0|
|θ1| = |-1/f 1|  |θ0|

Given the work in [2] Ch. 3, it's also easy to come up with simple transfer matrices for spherical or refraction interfaces, spherical mirrors, and thick lenses. We won't use those, but they're simple.

Express a paraxial system with elements described by 2x2 matrices, say M1, M2, M3, M4 in order from left to right. Then for ray r, M1*r is the ray after the first element, M2*M1*r, M3*M2*M1*r, and M4*M3*M2*M1 is the complete transformation and the system transfer matrix.

This ray-transfer technique is used to find other basic geometric properties describing an optical system, such as its nodal points and first and second principal planes. Not only that, but the four elements of the transfer matrix have discernible semantics, which can be illuminated by considering the physical meaning of setting each one to zero. For now, though, we're done.


Let's find some important distances for an example optical system. Consider three lenses of focal length 200, -50, and 50 mm. The first and second are separated by 100 mm. the second and third by 50 mm. The input plane is 200 mm. in front of the first lens. What is the (back, effective) focal length of the system? focal

We send a ray parallel to the axis into the 6-element system (translation, lens, translation, lens, translation, lens). The last ray we get is the one emerging from the last lens at X=350. It turns out to be (3, -.04)T. We're after the focal length g= y/tan(θ), (the base of the final triangle in the figure from X= 350 to X= 425). With the small angle assumption, the tangent of the angle is the angle, so this ray intercepts the axis at g= y/θ, or 75mm. We can then add a final, seventh, transfer of this distance, which should bring the ray down to the axis at the focal point, and that final system produced the figure above.

The axial image point of a system is where the image of a ray starting out at the origin crosses the optic axis. It's like the focal point, and calculated the same way, only using a different initial ray. It turns out that all the system's image points from an object at the origin will fall in the plane at that distance, so that is where an in-focus image will be formed.

For fun let's use a slightly different 3-lens system: lenses are spaced out by 300, 100, and 50 of focal lengths 200, -50, and 50. Shooting out a ray from the origin at angle .01 and calculating the axial image distance (it's 95 mm. out from the third lens) we can compose a final translation after the 3rd lens of 95 mm to get this plot: axial

With this same system, let's shoot three rays from some non-zero y height at three different angles: they should intersect at the same axial image plane as above. Nicely, all the rays from the same object point wind up at the same image point (coincidence?!) and we can use the y height at the image plane to get the ratio object-height/image-height, or linear magnification, as .2 (the image height is upside-down).


2. Ray-Casting and Systems of Linear Equations

The Problem and the Method

Ray-casting is the basic operation in rendering computer graphics representations into images. The idea sounds simple: given a point of view and direction of gaze (ray), compute the color and intensity of the light in that direction. In geometrical optics the wave nature of light is ignored, and to follow light-paths (rays) through optical systems we only need geometry, some algebra, and a few physical laws.

A basic case is: given a 3-D point of origin for the ray and its direction, where does it intersect a given plane in space? No different from the "line-intersect plane" problem, and in optics sometimes called "intersecting a ray with a plane mirror", because mirrors are the most common use of optical "flats".

As usual, represent points p, x, r 3-D by (x,y,z)T vectors. Let's represent directions α, β etc. in 3-D by (x,y,z)T unit vectors. They form a family that lives on a sphere of unit radius centered on the origin. You can see that every direction corresponds to a vector (point) whose head is on this Gaussian Sphere.

The standard way to intersect a vector with something is to stretch it out in its direction until it hits; the resulting length is all we need. So: Any point on a ray can be written
Ray: r = r0 + d α,
with d the length, r0 the ray's origin, α its direction. This vector equation represents three linear equations in x,y,z.

The plane equation is linear, and a 3-D version of the familiar line equation:
Plane: Ax + By + Cz + D = 0
Written like this, [A, B, C]T is a vector specifying a direction, in fact the direction normal to the plane. Scaling the whole equation so that this direction is a unit vector gives the new, scaled value of D a meaning: it's the perpendicular distance from the plane to the origin. We're going to simplify our lives right off by forcing our infinite mirror, or plane, to pass through the origin, so we can describe it with a linear equation: one strictly in (x,y,z) (D = 0: no pesky constant).
Plane through Origin: Ax + By + Cz = 0, or
(A/C)x + (B/C)y + z = 0.

Aside: Ray-casting actually wants to be described with projective geometry, a fascinating and pretty branch of mathematics. A (pinhole) camera projects the scene onto the sensor by raycasting, and we get the familiar perspective effects like the parallel railroad tracks seeming to come together out in the distance. This is a non-linear transformation: linear transformations preserve paralellism. Luckily, there is a neat trick for doing the needed calculations linearly (strictly matrix operations) only in a simple four-dimensional version of space (homogeneous coordinates). For a more general raycasting project, still pretty easy, see Pinholes and Beachballs. That project visualizes spheres since infinite planes are boring and non-infinite planar surfaces are very fiddly.

The equations Ray and Plane above are four linear equations. A solution to them gives the d at which the ray intersects the plane, and we're done. Preview: we write down the equations, do a quick massage, write them in matrix form, and the answer is obvious.

Ray origin: (x0,y0,z0).
Ray direction: 123).
3 Ray Equations, 1 Plane Equation:
x = x0 + α1d
y = y0 + α2d
z = z0 + α3d
(A/C)x + (B/C)y + z = 0.

We need four equations for unknowns x,y,z,d, but clearly we only really need three numbers, x,y,d say, from which we can easily solve for z using the Plane through Origin equation. Preemptively using that equation first, we get an expression for z to substitute into the third ray equation and we can rewrite the system as:
x0 = x - α1d
y0 = y - α2d
z0 = -(A/C)x - (B/C)y - α3d.

We could keep on substituting to solve the system, but let's express it as a matrix equation:
[x0,y0,z0]T = M [x,y,d]T , where M is
| 1             0   -α1 |
| 0             1   -α2 |
| -(A/C) -(B/C) α3 |

But hold it! If u = M v, then v = M-1 u, no? So we can solve for what we want (x,y,d) in terms of what we know or can easily compute ( x0,y0,z0, M-1). Then we ignore d, put x,y into the Plane through Origin equation, get z, and we've got our intersection-point.


  1. Tom Brown; Course Notes, Optics 211 2008, Part I.
  2. F.L. Pedrotti and L.S. Pedrotti; Introduction to Optics, Prentice-Hall 1987.
  3. Wikipedia.
  4. R. W. Ditchburn; Light, Wiley Interscience, N.Y. 1963.

Last Change: 05/02/2011: CB