Combining box filters, AA filters and diffraction: Do I need an AA filter?

I have been building up to this post for some time now, so it should not surprise you too much. What happens when we string together the various components in the image formation chain?

Specifically, what happens when we combine the square pixel aperture, the sensor OLPF (based on a 4-dot beam splitter) and the Airy function (representing diffraction)? First off, this is what the MTF curves of our contestants look like:

The solid black curve represents a combined sensor OLPF (4-dot beam splitter type) + pixel aperture + lens MTF (diffraction only) model. This was recently shown to be a good fit for the D40 and D7000 sensors. The dashed blue curve represents the MTF of the square pixel aperture (plus diffraction), i.e., a box filter as wide as the pixel. The dashed red curve illustrates what a Gaussian MTF (plus diffraction) would look like, fitted to have an MTF50 value that is comparable to the OLPF model. Lastly, the solid vertical grey line illustrates the remaining contrast at a frequency of 0.65 cycles per pixel, which is well above the Nyquist limit at 0.5 cycles per pixel (dashed vertical grey line).

Note how both the Gaussian and the OLPF model have low contrast values at 0.65 cycles per pixel (0.04 and 0.02, respectively), while the square pixel aperture + lens MTF, representing a sensor without an AA filter, still has a contrast value of 0.27. It is generally accepted that patterns at a contrast below 0.1 are not really visible in photos. That illustrates how the OLPF successfully attenuates the frequencies above Nyquist, but how does this look in a photo?

Ok, but how would it affect my photos visually?

I will now present some synthetic images to illustrate how much (or little) anti-aliasing we obtain at various apertures, both with and without an AA filter. The images will look like this:

The left panel is a stack of four sub-images (rows) separated by white horizontal bars. Each sub-image is simply a pattern of black-and-white bars, with both black and white bars being exactly 5 pixels wide (in this example). The four stacked sub-images differ only in phase, i.e., in each of the four rows the black-and-white pattern of bars is offset by a horizontal distance between 0 and 1 pixels in length.

The right panel is a 2x magnification of the left panel. Note that the third row in the stack is nice and crisp, containing almost pure black and pure white. The other rows have some grey values at the transition between the black and white bars, because the image has been rendered without any anti-aliasing.
These images are rendered by sampling each pixel at 2362369 sub-pixel positions, weighting each sampled point with the relevant point spread function.

The aliasing phenomenon known as frequency folding was illustrated in a previous post. When a scene contains patterns at a frequency exceeding the Nyquist limit (highest frequency representable in the final image), the patterns alias, i.e, the frequencies above Nyquist appear as patterns below the Nyquist limit, and are in fact indistinguishable from real image content at that frequency. Here is a relevant example, illustrating how a frequency of 0.65 cycles per pixel (cycle length of 1.538 pixels) aliases onto a frequency of 0.35 cycles per pixel (cycle length of 2.857 pixels) if no AA filter is present:

 

This set was generated at a simulated aperture of f/1.4, which does not attenuate the high frequencies much. Observe how the two images in the "No OLPF" column look virtually the same, except for a slight contrast difference; it is not possible to tell from the image whether the original scene contained a pattern at 1.538 pixels per cycle, or 2.857 pixels per cycle.

The "4-dot OLPF" column shows a clear difference between these two cases. If you look closely you will see some faint stripes in the 2x magnified version at 1.538 pixels per cycle, i.e., the OLPF did not completely suppress the pattern, but attenuated it strongly.

If we repeat the experiment at f/4, we obtain this image:

At f/4, we do not really see anything different compared to the f/1.4 images, except an overall decrease in contrast in all the panels.

Ok, rinse & repeat at f/8:

Now we can see the contrast in the "No OLPF" column, at 1.538 pixels per cycle, dropping noticeably. Diffraction is acting as a natural AA filter, effectively attenuating the frequencies above Nyquist.

Finally, at f/11 we see some strong attenuation in the sensor without the AA filter too:

You can still see some clear stripes (top left panel) in the 2x magnified view, but in the original size sub-panel the stripes are almost imperceptible.

Conclusion

So there you have it. A sensor without an AA filter can only really attain a significant increase in resolution at large apertures, where diffraction is not attenuating the contrast at higher frequencies too strongly. Think f/5.6 or larger apertures.

Unfortunately, this is exactly the aperture range in which aliasing is clearly visible, as shown above. In other words, if you have something like a D800E, you can avoid aliasing by stopping down to f/8 or smaller, but at those apertures your resolution will be closer to that of the D800. At apertures of f/5.6 and larger, you may experience aliasing, but you are also likely to have better sharpness than the D800.

Not an easy choice to make.

Personally, I would take the sensor with the AA filter.

Nikon D40 and D7000 AA filter MTF revisited

In an earlier post, I have shown some early MTF measurements taken with both the D40 and the D7000 at the centre of a Sigma 17-50 mm f/2.8 lens.
In this post, I revisit those measurements, presenting some new results for the D7000, and a model of the sensor AA filter (or OLPF).

Optical Low Pass Filters (OLPFs) in DSLRs

I have noticed previously that the measured MTF of the D40 was fairly close to a Gaussian, but that there were some small systematic discrepancies that were not accounted for. But how do you build an optical filter that has a Gaussian transfer function?

It turns out that the OLPFs in DSLRs are not really Gaussian, but something much simpler: beam splitters. A slice of Lithium Niobate crystal, cut at a specific angle, presents a different index of refraction depending on the polarization of the incident light. Let us assume that a beam of horizontally polarized light passes through straight, for the sake of the argument. Vertically polarized light, on the other hand, refracts (bends) as it enters the crystal, effectively taking a longer path through the crystal. As the vertically polarized light leaves the crystal, it refracts again to form a beam parallel to the horizontally polarized beam, but displaced sideways by a distance dependent on the thickness of the crystal.

Using a single layer of Lithium Niobate crystal, you can split a single beam into two parallel beams separated by a distance d, which is typically chosen to match the pixel pitch of the sensor. Since this happens for all beams, the image leaving the OLPF is the sum (average) of the incoming image and a shifted version of itself, translated by exactly one pixel pitch.

If you stack two of of these filters, with the second one rotated through 90 degrees, you effectively split a beam into four, forming a square with sides equal to the pixel pitch (but often slightly less than the pitch, to improve resolution). A circular polariser is usually inserted between the two Lithium Niobate layers to "reset" the polarisation of the light before it enters the second Niobate layer.

Combining pixel aperture MTF with the OLPF MTF

So how does this beam-splitter effect the desired blurring? We can compute the combined PSF by convolving the beam-splitter impulse response with the pixel aperture impulse response (a box function).

The beam splitter is represented as four impulses, i.e., infinitely thin points. Basic Fourier transform theory tells us that the Fourier transform of an impulse is just a cosine function, so the MTF of the beam splitter will be the sum of four cosines.

The "default" 4-way beam splitter impulse response filter can be denoted as the sum of four diagonally-placed impulses, i.e.,

f(x,y) = δ(x-d, y-d) + δ(x+d, y+d) + δ(x-d, y+d) + δ(x+d, y-d)

where δ(x,y) represents a 2D Dirac delta function, which is non-zero only if both x and y are zero, and d represents the OLPF split distance. More sophisticated OLPF designs are possible (e.g., 8-way splitters), but the 4-way designs appear to be popular. In my notation here, the distance between two beams would be 2d; this is to accommodate my usual notation of a pixel being defined over the area [-0.5, -0.5] to [0.5, 0.5].


The degree of blur is controlled by the d parameter, with d = 0 yielding no blurring, and d = 0.5 giving us a two-pixel-wide blur. Since d can be varied by controlling the thickness of the Lithium Niobate layers, a manufacturer can fine-tune the strength of the OLPF for a given sensor.

Convolving a square pixel aperture with four impulse functions is straightforward: just sum four copies of the box filter PSF, each shifted by d in the required direction. For d = 0.375, we obtain the following PSF:

in two dimensions, or

in three dimensions. Not exactly a smooth PSF, but then neither is the square pixel aperture.

Simulating the combined effects of diffraction, OLPF and pixel aperture

We can derive the combined PSF directly by convolving the diffraction, OLPF and pixel aperture PSFs. Note that this combined PSF is parameterized by wavelength, lens aperture, pixel pitch, and OLPF split distance. For a wavelength of 0.55 micron, an aperture of f/4, a pixel pitch of 4.73 micron and an OLPF split distance of 0.375 pixels (=2.696 micron), we obtain the following PSF:

in two dimensions, or

in three dimensions. Notice how diffraction has smoothed out the combined OLPF and pixel aperture PSF (from previous section).

Using this PSF we can generate synthetic images of a step edge. The MTF of this synthetic edge image can then be compared to a real image captured with a given sensor to see how well this model holds up.

Enter the Death Star

You can certainly measure lens or sensor sharpness (MTF) using any old chart that you printed on office paper, but if you are after resolution records, you will have to take a more rigorous approach. One of the simplest ways of obtaining a reasonably good slanted edge target is to blacken a razor blade with candle soot. This gives you a very straight edge, and very good contrast, since the soot does not reflect much light.

That leaves only two other questions: what do you use as a background against which you will capture the razor blade, and how do you illuminate your target?

Previously I have used an SB600 speedlight to illuminate my razor blade, which was mounted on a high-grade white paper backing. This produced reasonably good results, but it did not maximize contrast because the flash was lighting the scene from the front. There is also a possibility that the D7000 cycles its mirror when using a flash in live view mode, which could lead to mirror slap. So the flash had to go.

In its place I used a home made integrating sphere, which I call the "Death Star" (sorry George, please do not sue me). Here is what it looks like with the razor blade mounted over the integrating sphere's port:

Yes, I know the highlight is blown out ....

Why use an integrating sphere? Well, the integrating sphere produces perfectly uniform diffuse lighting, which is ideal for creating a uniform white backdrop for the razor blade.

In addition, my home made integrating sphere produces enough light to require a shutter speed of 1/200 s at f/4, ISO100 to prevent blowing out the highlights. This is perfect for minimizing the influence of vibration (although higher shutter speeds would be even better).

Using this set-up, I then capture a number of shots focused in live view, using a remote shutter release. It appears that with this target, the AF is really accurate, since I could not really do much better with manual focus bracketing. One day I will get a focusing rail, which will make focus bracketing much simpler.

Lastly, all MTF measurements on the razor edge were performed using MTF Mapper. The "--bayer green" option was used to measure MTF using only the green photosites, thus avoiding any problems with Bayer demosaicing. The raw files were converted using dcraw's "-D" option.

D40 MTF and OLPF model

Here is the MTF plot of a D40 razor image captured at f/4 (manually focus bracketed), Bayer green channel only:

D40 green channel MTF, Sigma 17-50 mm f/2.8 at f/4

The D40's PSF was modelled by convolving the square pixel aperture (7.8 micron wide), the 4-point beam splitter (d=0.375 pixels), and the Airy function (f/4). A synthetic image of a step edge with this PSF was generated, and measured using MTF Mapper.

Purely judging the match between the model and the measured MTF by eye, one would have to conclude that the model captures the interesting parts of the MTF rather well. The measured MTF is slightly below the model, which is most likely caused by a smidgen of defocus.

The fact that the model fits so well could also be taken to imply that the Sigma 17-50 mm f/2.8 lens is relatively aberration-free at f/4, i.e., it is diffraction limited in the centre.

MTF50 resolution came in at  0.334 cycles per pixel, or 42.84 line pairs per mm, or 671 line pairs per picture height.

D7000 MTF and OLPF model

Here is the MTF plot of a D7000 razor image captured at f/4 (AF in live view), Bayer green channel only:

D7000 green channel MTF, Sigma 17-50 mm f/2.8 at f/4

The D7000's PSF was modelled by convolving the square pixel aperture, the 4-point beam splitter (d=0.375 pixels), and the Airy function (f/4). A synthetic image of a step edge with this PSF was generated, and measured using MTF Mapper.

The measured MTF does not fit the model MTF quite as well as it did in the D40's case. Given that the physical linear resolution is 60% higher, it is correspondingly harder to obtain optimal focus. The shape of the measured MTF relative to the model MTF is consistent with defocus blur.


The actual resolution figures are impressive: MTF50 is 0.304 cycles per pixel, or about 64 lp/mm, or equivalently, 992 line pairs per picture height.

If the model is indeed accurate, it would mean that the D7000 can theoretically obtain resolution figures close to 68 lp/mm at f/4 in the green channel, provided that the lens is purely diffraction limited, and perfectly focused.

Summary

Perhaps this is not such a surprising result, but it appears that Nikon is using the same relative strength of AA filter in both the D40 and the D7000; this can be deduced from the fact that both the D40 and the D7000 OLPF models fitted best with an OLPF split distance of 0.375 pixels.


The somewhat unexpected result, for me at least, was that the MTF shape is so sensitive to perfect focus. Specifically, it seems that the first zero of the MTF curve, at around 0.6875 cycles per pixel, is not readily visible unless focus is perfect. The zero is quite clear in the D40 curve, but not quite so visible in the D7000 curve. You are extremely unlikely to achieve this kind of focus in real world photography, though.

References

1. Russ Palum, Optical Antialiasing filters, in Single-Sensor Imaging: Methods and Applications for Digital Cameras, Edited by Rastislav Lukac, CRC Press 2008

D800E versus diffraction

In a previous post (here), I have illustrated how diffraction through a circular aperture can be modelled either in the spatial domain as a point spread function (PSF), or in the frequency domain as a modulation transfer function (MTF). I will now put these models to use to investigate the influence of diffraction on the resolution that can be achieved with the D800E at various apertures.

Simulating the effect of diffraction

I will not go into the maths behind the diffraction MTF; this was discussed in another post (here). For now, it is sufficient to understand that we can combine the diffraction MTF with the sensor's MTF through multiplication in the frequency domain.

Assume for the moment that the D800E effectively does not have an AA filter (in practice, this might not be entirely true, i.e., the D800E may just have a very weak AA filter compared to other cameras). This allows us to model the pixel's MTF curve as a sinc(x), as was shown in a previous post. Next, we assume that the lens is diffraction limited, i.e., the other lens aberrations are negligible, and thus the lens MTF is just the diffraction MTF.

For a D800(E) pixel pitch of 4.88 micron, and an aperture of f/8, we obtain the following combined MTF curve:

D800E combined MTF at f/8

The dashed grey curve represents the sensor's MTF, and the black curve represents the diffraction MTF. The blue curve is the product of these two curves, and represents the combined diffraction-and-sensor MTF.

At f/8, our peak MTF50 value will be 0.344 c/p, or 70.4 lp/mm. Note that this is still higher than what I measured on a D7000 at f/5, which peaked at about 0.29 c/p (61 lp/mm), but the D7000  has an AA filter. 

Moving to even smaller apertures will cost us resolution, thus at f/11 the curve looks like this:

D800E combined MTF at f/11

At f/11, MTF50 peaks at only 0.278 c/p, or 57 lp/mm. This is still extremely crisp, although you might barely be able to see the difference compared to f/8 under ideal conditions. Pushing through to f/16:

D800E combined MTF at f/16

Note how close the combined MTF curve and the diffraction MTF curve have now become; this indicates that diffraction is starting to dominate the MTF curve, and thus also resolution. At f/16, MTF50 has dropped to 0.207, or about 42.3 lp/mm, which is not bad, but quite far from the 70 lp/mm we achieved at f/8.

What about going in the other direction? Here is what happens at f/5.6:

D800E combined MTF at f/5.6

MTF50 now reaches 0.412 c/p, or 84.4 lp/mm. At f/4 (not shown as a plot) we get 0.465 c/p (95.3 lp/mm), and so on. Below f/4 we will start seeing the residual aberrations of the lens take over, which will reduce effective resolution. I have no model for those yet, so I will stop here for now.

Ok, so I will go one step further. Here is the MTF plot at f/1.4, but keep in mind that for a real lens, other lens aberrations will alter the lens MTF so that it is no longer diffraction limited. But this is what it would have looked like if those aberrations were absent:

D800E combined MTF at f/1.4

Off the charts! MTF50 will sit at 0.557 c/p, or 114.1 lp/mm. The pixel MTF and the combined MTF are now very similar, which is to be expected, since diffraction effects are now almost negligible. Now if only they could build this lens ...

In closing

These results seem to support the suggestions floating around on the web that the D800E will start to visibly lose sharpness after f/8, compared to what it achieves at f/5.6. But this does not mean that f/11 is not sharp, since 57 lp/mm is not something to be sneezed at! Even more importantly, there is no "magical f-stop" after which diffraction causes the resolution to drop; diffraction will lower resolution at all f-stop values. The balance between diffraction blur and blur caused by other lens aberrations tends to cause lens resolution to peak at a certain aperture (around f/4 to f/5.6 for many lenses), but even at f/1.4 you will lose resolution to diffraction, just not a lot.

There are also some claims that the D700 was usable at f/16, but now suddenly the D800E will not be usable at f/16 any more. This is not true. If we compare a hypothetical D700E with our hypothetical D800E above, we see that the D800E attains an MTF50 value of 42.3 lp/mm at f/16, and the hypothetical D700E would reach only 37.2 lp/mm.

The real D700 has an AA filter. If we approximate the strength of this filter as a Gaussian with a standard deviation of 0.6246, then the D700 would only reach an MTF50 of 25.6 lp/mm at f/16. A similar approximation of the AA filter for the D800 would produce an MTF50 of 34.4 lp/mm at f/16. So the D800 (or D800E) will always capture more detail than the D700 at all apertures. The D800E is perfectly usable at f/16, and more so than the D700. 

[Incidentally, the diffraction + Gaussian AA filter approximation used here appears to be quite accurate. Roger Cicala's Imatest results on the D800 and D700 with the Zeiss 25 mm f/2 (see here) agree with my figures. From Roger's charts, we see the D800 at f/5.6 achieves 1200 lp/ph, or about 50.06 lp/mm, compared to my figure of 50.7 lp/mm. The D700 at f/5.6 attains roughly 750 lp/ph (31.38 lp/mm) in Roger's test, and my model predicts 31.9 lp/mm.]

The catch, though, is that the D700's MTF50 at f/16 is 0.216 c/p (25.6 lp/mm), whereas the D800's MTF50 at f/16 is 0.168 c/p (34.4 lp/mm). The apparent per-pixel sharpness of the D700 will exceed that of the D800 at 100% magnification on-screen. If you view them at the same size, though, the D800 will be somewhat sharper.

An introduction to diffraction

Incorporating diffraction

A while ago I posted an article on the MTF of a sensor without an AA filter in the absence of diffraction (here if you've missed it). By not including the effects of diffraction, we could see how the lack of an AA filter produced significant aliasing, which shows up visually as false detail.

Now it is time to add the effects of diffraction back in. I do not have the required background in physics to explain this from first principles, but we can probably blame diffraction on the wave properties of light. If we take a beam of light and pass it through a narrow vertical slit, and then project it onto a surface, we do not see quite what we expected.

Instead of seeing this:

we actually see this

Intensity pattern of diffraction through narrow slit

Of course, to illustrate this principle I had to compress the intensity of the diffraction image quite a bit, since the secondary stripes are much dimmer than the main one. Making the slit narrower causes the intensity pattern to spread out more.
The intensity pattern that we observe is in fact sinc²(x), and looks like this:

Intensity profile of narrow slit, sinc²(x)

Recall that we first observed the sinc(x) function, which is defined as

sinc(x) = sin(x) / (x),

when we computed the Fourier transform of the box function. Note that a cross-section of the vertical slit (along constant y values) is in fact a box function. Why is the intensity function sinc²(x) rather than sinc(x) ? Again, I lack the physics background to explain this properly, but it is squared because we are dealing with incoherent light, which is a safe assumption for typical photographic light sources. Coherent light sources, such as lasers, will have a sinc(x) diffraction pattern upon passing through a narrow slit.

If we have a square hole (aperture) rather than a vertical slit, we will observe this intensity pattern:

Intensity pattern of square aperture

which is simply sinc²(x)*sinc²(y).

If we want to introduce a circular aperture, like the one in a real lens, then we can no longer treat x and y separately. A circular aperture produces an intensity pattern that looks like this:

Intensity pattern of circular aperture (Airy pattern)

This pattern is also known as the Airy pattern, with the bright central circular area (before first black ring) known as the Airy disc. Again, a smaller aperture causes the size of this disc to increase, and the entire pattern to spread out.

Although x and y can no longer be treated independently, we can transform them to polar coordinates to obtain r = sqrt(x² + y²), which yields the intensity function jinc²(r) where

jinc(r) = 2*BesselJ(r)/(r)

and BesselJ is a Bessel function of the first kind, of order 1 (see Wikipedia). Since this function is radially symmetrical (as can be seen in the intensity image above) we can plot it simply as a function of r: 

Intensity profile (radial) of circular aperture, jinc²(r)

The result is that the circular aperture behaves somewhat like a square aperture, except that instead of observing a sinc²(x) function, we instead observe a jinc²(r). Note that the side-lobes of the jinc²(r) function are much smaller than those of the sinc²(x) --- you can probably guess what effect this will have on MTF.

Diffraction MTF

So what does the MTF curve of a diffraction pattern look like? For the first case, namely the MTF of the diffraction pattern of a vertical slit, we can perform the analysis in one dimension, to obtain the following curve:

MTF of sinc²(x) PSF

If you have guessed that it looks an awful lot like a straight line, then you guessed correctly. Recall that the intensity function of the vertical slit is sinc²(x). The Fourier transform of sinc(x) is a box function (see here for a recap), and vice versa, so the MTF of sinc(x) is just a box function in the frequency domain. Since multiplication in the time (or spatial) domain is equivalent to convolution in the frequency domain, we can see that the Fourier transform of sinc²(x) must be the convolution of a box function with itself in the frequency domain. Lastly, recall that the convolution of a box filter with itself (once only) yields a triangle function, and we are only plotting the positive frequencies in the MTF plot, i.e., you can mirror the MTF around the y axis to see the triangle more clearly, like this:

MTF of sinc²(x) PSF, showing negative and positive frequencies

This result extends to the MTF of the diffraction pattern of a square aperture, since x and y can simply be treated independently. As you have probably guessed, the 2D MTF of the square aperture looks like this:

MTF (left) and PSF (right) of square aperture

Note that this MTF shape can also be obtained through the convolution of two 2D box functions.

You may be tempted to treat the circular aperture in the same way, but trying to compute the MTF in one dimension will not give you the correct result. If you can visualise the MTF of the square aperture as the convolution of a box function with itself, then you can visualise the MTF of the circular aperture as the convolution of two "pillbox" functions that look like this:

A pillbox function

The convolution of two of these pillbox functions will give you the MTF on the left:

MTF (left) and PSF (right) of circular aperture

which is called the Chinese hat function, or chat(x). We can take a cross-section through the centre of this 2D function to obtain this plot:

chat(f) = MTF of jinc²(r)

Note that this function is very nearly a triangle function --- except that it curves a bit. This should not be totally surprising, since you can imagine that you could fit a square box function inscribed in the cylinder formed by the pillbox function.

If we focus only on the positive half of the frequency spectrum, the function can be expressed as

chat(s) = 2/π * (acos(s) - s*√(1 - s²))

which is the form that you will usually see in optics literature. Now we know what the MTF curve looks like, so we can obtain the corresponding point spread function by computing the inverse Fourier transform of the MTF curve.


If you perform this inversion in two dimensions, and then take a cross-section through the result, you will obtain the known intensity function, jinc²(r).
Right. Now sample the jinc²(r) function in one dimension, and apply the Fast Fourier Transform (FFT) to obtain the 1D MTF:

Trying to obtain chat(f) by computing FFT of jinc²(r)

Whoops! That does not work as expected. You have to perform the calculations in 2D for this to work. It only took me a (tortuous) weekend to discover this useful fact. If only someone had warned me ...

Other useful facts regarding diffraction

Just as with a Gaussian, the width of the MTF function is inversely proportional to the width of the PSF, i.e., smaller apertures produce larger diffraction patterns.

Diffraction is wavelength-dependent. If jinc²(r) is your PSF, then r should be expressed as

r = π * d / (λ * N)

where d is your radial distance from the optical axis, as measured in the focal plane. If you are modelling the diffraction PSF in terms of pixels, then you can express d in pixels, but you must adjust λ accordingly, such that

r_p = π * d_p / ( (λ/Δ) * N)

where Δ denotes your pixel pitch, expressed in the same units as λ. For example, the Nikon D800 has a pixel pitch of 4.88 micron, and green light has a wavelength of 550 nm, or 0.55 micron. This means that the diffraction PSF expressed in D800 pixel-units would be

r_p ≈ π * d_p / ( 0.112705 * N).

If we want to obtain the equivalent function in the frequency domain, i.e., the diffraction MTF, we use the scale

s_p = (λ/Δ) * N * f_p

where f_p  is the frequency scale in cycles per pixel, and s_p is normalized frequency we plug into the equation

chat(s_p) = 2/π * (acos(s_p) - s_p*√(1 - s_p²)).

For the D800 pixel size, we would use

s_p ≈ 0.112705 * N * f_p

In a future post, some examples of MTF curves, incorporating diffraction, will be investigated using the D800E as example.

Steve, Otto and the green-screen

This post is not about MTF Mapper or lens sharpness. This story is about Steve.

And about chroma keying, commonly called "green screen" or "blue screen" effects. This technique is mostly needed by Hollywood. You cannot have an actor drive a car while having a spirited argument with another actor, all the time whilst driving in actual traffic. Nor can you quickly pop out to outer space to get a few shots of the starship Enterprise.

The typical solution is to capture the footage of your foreground objects, typically the actors or your model starship, against a uniformly coloured green background. Then you assign a transparency value (called an alpha channel) to each pixel in your image so that where the green background is showing through, you have full transparency. Hopefully the pixels that comprise the actor will have zero transparency, but you have to allow for partial transparency in the actor's hair, for example.

The problem can thus be stated as follows:
The green screen problem: Given a shot of the foreground object (actor or whatnot) against a uniform green (or blue) background, compute the correct transparency value for each pixel in the image.

Of course, you have to ensure that your foreground object does not have any parts that are close to the background colour, or these parts will become transparent.

James Blinn (and co-author A.R. Smith, in a paper titled "Blue screen matting") have shown that the green screen problem is in fact an underdetermined problem, meaning that you have more variables to determine than you measurements. In the simplest case, you know the background colour, in RGBA format, (R_b, G_b, B_b, 1.0). Note that the background is fully opaque, so its alpha is exactly 1.0. The foreground object's colour is (R_f, G_f, B_f, A_f), where it is assumed that the RGBA values have already been pre-multiplied by the alpha (opacity) value. This means that the colour of a pixel in our captured image is simply

(R_i, G_i, B_i, 0.0) = (R_f, G_f, B_f, A_f) + (1-A_f)*(R_b, G_b, B_b, 1.0)
We can break this down to give

R_{i =} R_{f +} (1-A_f)*R_b
and similarly for the other two colours. This gives us three equations with four unknown values (R_f, G_f, B_f, A_f), which cannot be solved uniquely.

So, despite the fact that the maths shows you that green-screen techniques cannot work, Hollywood insists on using it! Of course, if you are willing to add some assumptions (constraints) to the above equations, you can find some unique solutions. One of these constraints is that the foreground colour is distinct from the background colour, e.g., no green actors. Even under such limiting assumptions we often find that the results are not really all that great; just think of the flying-on-the-broomstick scene from Harry Potter and the Order of the Phoenix movie.

Blinn and Smith have shown that unique solutions can be obtained by capturing the foreground object against two different backgrounds. You also have to capture an image of each background without any foreground objects. Although this is somewhat of a bother, it allows you to obtain much better results in practice. Remember that the foreground object opacity, A_f, is same in the two images against different backgrounds, but that the apparent foreground object colour (R_f, G_f, B_f) will not necessarily be the same if the object is not fully opaque. Combining the equations from both images thus gives us six equations with six unknowns, which fortunately will have a unique solution, provided that the two backgrounds really are different for each pixel.

The solution suggested by Blinn and Smith goes as follows:
Let D = (R_i,1 - R_i,2, G_i,1 - G_i,2, B_i,1 - B_i,2)
where  R_i,1 denotes the red value of our pixel of interest in image 1, and R_i,2 denotes the red value of the same pixel in image 2. In other words, D is the difference between the two input images containing the foreground objects against the two different backgrounds.
Similarly, let E = (R_b,1 - R_b,2, G_b,1 - G_b,2, B_b,1 - B_b,2)
denote the difference between the two background-only images, without the foreground objects, i.e., the pair of background-only shots you have to capture as part of this method.
The alpha value for a pixel is thus given by
A_f = 1 - (D · E) / (E · E)
where  (D · E) denotes the dot product between the two vectors.

And that is all there is to this method; you can reconstruct the foreground object's red component as
R_f = (R_i,1 - (1.0 - A_f)*R_b,1)
and similarly for the other two colours.
Now you can recombine the foreground objects with a new background as
R_n = R_f + (1.0 - A_f)*R_k
where R_k denotes the red component of the new background image. Just apply the same pattern for blue and green components, and you are done. Remember, R_f is pre-multiplied with A_f, which explains why the re-blending with a new background appears in this form.

How well does this work? I called in the help of Steve and Otto to demonstrate. Here is the first shot against background 1 (click for a larger version; applies to all images in this post):

Input image 1, shot against background 1

This image contains several interesting things. Note that there is quite a bit of green in the foreground objects (duck's head, for example). The champagne flute is also a bit tricky, because we can see most of the background right through it. Otto's fur would have presented endless trouble if you tried to remove the background manually in Gimp, for example. Yes, the dog is called Otto (Steve is the lime).

And the second shot:

Input image 2, shot against background 2

Lastly, a shot of each of the backgrounds, after physically removing the foreground objects from the scene.

background 1

background 2

I pushed these four images through a quick-and-dirty implementation of Blinn and Smith's method that I cobbled together in C++ using OpenCV. Here is the resulting alpha mask produced by the program without any manual intervention:

Alpha mask produced using Blinn and Smith's method

Note how solid Steve's alpha mask is --- the fact the he's green, and one of the backgrounds happened to be green, made no difference whatsoever. Otto's interior is also completely opaque, but on the edges we see some fine details like these:

100% crop of the alpha mask near Otto's head

Notice how the fine fibres are partially transparent (top right), and how the interior gaps in the fur is also partially transparent (bottom left-ish).

The champagne flute has also been extracted quite nicely.

Now for the final test: to recombine the foreground objects with a new background. Here is the background image:

PovRay's benchmark scene should do nicely as a new background

And here is the blended result:

Recombining the extracted foreground with a new background

The new composition is essentially flawless. There are no tell-tale fringes or other signs that are commonly seen with green-screen techniques.

Here is a close-up of Otto's head after blending:

100% crop of Otto's head after blending

Of course, if you concentrate a bit, you will see something is amiss. Although the champagne flute has been composited perfectly, it is quite clear that Steve (the lime) is correctly refracted through the right-hand side of the glass, but that the PovRay logo in the background suffered no distortion on the left side of the glass. This is obviously a shortcoming of all such recomposition techniques, so I guess the only solution is avoid strongly-refractive foreground objects.

Limitations

Obviously it is somewhat inconvenient to have to shoot the subject against two different backgrounds. It will be downright impossible to use this technique with a toddler, for example. An adult might be able to hold still enough, provided you use something like a television to display your background like I did above.

This technique is fine for inanimate objects, though.

Oh yes, of course you have to keep the camera very, very still during the whole capture process. I used a sturdy tripod + mirror lock-up + IR shutter release. In theory, you could use image-to-image registration to correct for camera movement, but I have not tried it yet for this problem.

I recently saw another paper (published in the proceedings of SIGGRAPH'06) by McGuire and Matusik where they extend this technique to live video. The secret is to use a polarised background material, which will obviously present different images through different polarisations of the light. A special camera uses a prism to send the incoming light to two different sensors according to its polarisation. Although this means the technique falls into the category of "requires special tools", it is still pretty cool. And it produces perfect results, just like the method above. So why didn't they use this for the dreaded broomstick scene in Harry Potter and the Order of the Phoenix?

Notes:

No animals were harmed during the production of this article. Steve, however, did not make it. That was some good lemonade, though!

Pixels, AA filters, Box filters and MTF

What happens to a sensor's MTF curve when you remove the Anti-Aliasing (AA) filter (also called the Optical Low-Pass Filter, or OLPF)? I have chosen a rather long path to reach the answer, so this article is a bit longer than ideal, but I hope that you will agree that it is worth the effort to slog through it :)

To simplify the discussion, I am going to limit myself to grayscale sensors, e.g., the Leica M-Monochrom, or a lab camera such as a Prosilica. With no Bayer filter, there is little reason to add an AA filter to a sensor ... or is there?

First up we have to consider how such a grayscale sensor works. The sensor consists of a grid of photosites. The distance between the centres of two neigbouring photosites (left-right or top-bottom neighbours) is called the pixel pitch. In the ideal case, the pixel can be thought of as a little square with side lengths equal to the pixel pitch. In the real world, there has to be some space between adjacent pixels, and some space for circuitry to facilitate the read-out of a photosite. The fraction of the pixel that is actually collecting photons is called the fill factor. In the old days, this meant that only about 50% or so of your total sensor area was actually collecting photons, but with advances such as back-side illumination and microlens arrays we are doing much better today at about 90% fill factor.

Real world notwithstanding, it is convenient to think of a photosite as a little square, and a sensor as a grid of tightly packed photosite with practically no gaps between them.

This abstraction yields a very specific point spread function shaped like this:

The intensity of a grayscale sensor pixel is thus the sum of all the light collected over the surface of the photosite.

Enter the box filter

A 1D box filter is a filter that looks like this:

Clearly, the point spread function of an ideal square photosite is just the 2D version of a box filter.

The box filter pops up many times in signal processing contexts. It also crops up implicitly in unexpected places, but more on that later.

At this point I have to assume that you are familiar with convolution (check out the Wikipedia entry for a refresher). The convolution of a box filter f with itself (f * f) yields a triangular function like the black curve in this plot:

The red curve represents repeated convolution of a box filter with itself three times (((f * f) * f) * f), and the green curve 5 times. And yes, you guessed correctly, it does look a little like a Gaussian distribution. This should not surprise you any more.

I shall state without proof that repeated convolution of a box filter of width 1 pixel with itself does indeed converge on a Gaussian function with standard deviation of 1/√2 ≈ 0.7071.

The unexpected upshot (well, I certainly was not expecting it the first time) of this is that you can approximate a Gaussian blur by repeatedly blurring an image with a box filter. How do you blur an image with a box filter?

Well, convolution of an image with a 2D box filter that is 3x3 pixels in size is equivalent to replacing each pixel's value with the unweighted average of the group of 3x3 pixels centred at the pixel you are replacing. In this case, repeated convolution with a 3x3 2D box filter will approximate a Gaussian blur with a standard deviation of 3/√2 ≈ 2.12.

This observation, namely that repeated convolution with a box filter converges on a Gaussian, will be recycled in a later post, so do not page it out of your main memory just yet ...

Area-weighted Anti-Aliasing in computer graphics rendering

Rendering an image of a black square on a white background is relatively straightforward. You scan through each pixel, and test whether the current pixel is inside the square (pixel set to black) or outside (pixel set to white). This produces perfect images as long as the square is aligned with pixel boundaries, which is rather limiting.

If you rotate your black square relative to the pixel grid, a third possibility arises: a pixel might fall partially inside the square, and partially outside. The most intuitive solution to the problem is to simply set the intensity of the pixel proportional to the degree to which it is inside the square. In computer graphics terms, you clip the square to the boundaries of your pixel, and measure the area of the clipped polygon. Here is a picture:

This allows us to draw the black square with much smoother-appearing edges, since we are effectively using our radiometric resolution (different shades of gray) to compensate for our inadequate spatial resolution (large pixels) --- a practice commonly referred to as Anti-Aliasing in computer graphics.

We can also view this from a sampling perspective. First, we generate a number of (x,y) coordinates within each pixel, preferably well-spread throughout the current pixel. By counting how many of these coordinates are inside the square too, we have estimated the fraction of the current pixel that is covered by the square, which, in the limit, is identical to the fragment area that we obtained by clipping the square to the pixel boundaries above.

Lastly, we can consider the case where we sample at coordinates outside of the current pixel too, but that we apply a weight to each sample. If the weight is zero everywhere outside the current pixel, and exactly 1.0 inside the current pixel, then we can say that our weights represent the point spread function of the current pixel, which happens to be a 2D box function by construction.

This shows us that there are three equivalent methods of obtaining the correct intensity of a pixel that is only partially covered by the square:

1. We can compute the area of intersection between the current pixel's boundaries and the square, or
2. We can perform unweighted sampling at points spread uniformly throughout the current pixel, or
3. We can perform weighted sampling at arbitrary points, and weight each point relative to the box filter centred at the current pixel (obviously we will improve our efficiency by only sampling close to the current pixel).

By now you can see the link between Area-weighted anti-aliased rendering and an image sensor with square pixels. If our sensor is observing a knife edge target (or just a black rectangle on a white background), then the intensity of each photosite will be proportional to the area covered by the black square, i.e., the sensor is implicitly applying a box filter while sampling the real-world scene.

MTF of a box filter

The point spread function of a non-Bayer sensor without an AA filter is simply a box filter of width equal to the pixel pitch. The point spread function of a synthetic image generated with an area-weighted sampling algorithm is also a box of width equal to the pixel size, which means we can use the one to study the other. So what does the MTF of this box function look like?

Again, I shall state without proof that the Fourier transform of the box function is sinc(f), or sin(f)/f. Thus, if our point spread function is a box function, then our MTF will simply be sinc(pi*f)/(pi*f), which looks like this:

A couple of things are important about this MTF curve (black curve):

• Recall that in this case, a frequency of 0.5 cycles per pixel implies that our image contains a pattern of alternating black-and-white stripes that are exactly one pixel wide; together, one black stripe and one adjacent white stripe makes one cycle. This is the highest spatial frequency that can be represented correctly in our image --- if you try to make the stripes less than one pixel wide, then clearly you will not be able to preserve details exactly any more.
• Note that the contrast drops to zero at twice the Nyquist frequency (1 cycle per pixel). When our stripes are exactly half a pixel wide, then the black and white pattern cancels exactly, leaving the entire image at 50% grey, hence zero contrast.
• Also note that the box filter MTF curve is nonzero between 0.5 cycles per pixel and 1 cycle per pixel, which is highly undesirable. Aliasing is when a high frequency component masquerades as a lower frequency component (more on this below). As mentioned above, if the width of each stripe in a black-and-white pair falls below 1 pixel, then we can no longer represent this accurately in our image. 
• At each integer k cycles per pixel, we can fit k pairs of black-and-white stripes inside one pixel, again cancelling each other and producing a 50% grey image. This is exactly where the MTF curve drops to zero each time.
• The dashed grey curve is the MTF of a Gaussian PSF with a standard deviation 0.568 pixels (MTF50=0.33), included for comparison.

The MTF curve of the box filter describes exactly the MTF of a synthetic image rendered using Area-weighted AA (as described above), and it will be a good model of the MTF of a grayscale image sensor in the absence of the lens MTF (i.e., no diffraction).

Aliasing in practice

I have alluded to the fact that bad things happen when your scene contains detail that occurs at a higher frequency that what the sensor can capture, i.e., when you have details smaller than the pixels of your sensor. The type of aliasing I would like to highlight here is folding, which is when the high frequency information wraps around and re-appears as lower frequency information. To understand what it looks like, consider the following image:

This image does not exhibit any aliasing; it is merely to establish what we are looking for here. Firstly, the left panel is a stack of four sub-images (rows) separated by white horizontal bars. Each sub-image is simply a pattern of black-and-white bars, with both black and white bars being exactly 5 pixels wide. The four stacked sub-images differ only in phase, i.e., in each of the four rows the black-and-white pattern of bars is offset by a horizontal distance between 0 and 1 pixels in length.

The right panel is a 2x magnification of the left panel Note that the third row in the stack is nice and crisp, containing only pure black and pure white. The other rows have some grey values at the transition between the black and white bars, because the image has been rendered with box-filtered anti-aliasing.

Here is the same image again, repeated to simplify comparisons further down:

Box filtered, bars are 5 pixels wide

Note that the edges of the bars are exhibiting the classical box-filtered anti-aliasing patterns; depending on the exact position of the bars relative to the pixels (i.e., the differences between the four rows), we see different shades of grey at the transition, between the four rows. Contrast that with the same pattern, but rendered using a Gaussian filter (the one depicted in the dashed grey curve in the MTF plot above), rather than a box filter, to reduce aliasing:

Gaussian filtered (standard dev. = 0.568 pixels), bars are 5 pixels wide

Here we can see that the edges of the bars are noticeably more blurry, but they do appear much smoother than the version with box filtering above. Also note that the four rows now look much more alike, which is an improvement on the box filtered version.

Now for the promised frequency folding. If we define a cycle as one black bar followed by one white bar, then the 5-pixel wide bars give us a cycle, or period, of 10 pixels. Frequency is simply 1 over period, or 1/10 cycles per pixel. Since the frequency is less than one cycle per pixel, we know that we can accurately represent the bars at this frequency. We know that 0.5 cycle per pixel corresponds to the Nyquist limit, i.e., the highest frequency that can be represented, which corresponds to a pattern of a 1-pixel wide black bar followed by a 1-pixel wide white bar. Frequency folding dictates that a pattern with a frequency of 1 - (1/10) = 0.9 cycles per pixel will be aliased with a frequency of 1/10 cycles per pixel. We know that 0.9 cycles per pixel corresponds to a cycle of 1.1111' pixels in length, which implies that the bars will have a width of 0.5555' pixels, which we also know cannot be represented accurately, since the bars are smaller than the pixels. If we render the same bar pattern with box filtering, but choosing a bar pattern in which the bars are 0.5555' pixels wide, we get this image:

Box filtered, bars are 0.5555' pixels wide

Not what you were expecting? One would have expected that the bars disappear into a uniform grey patch, since we know that 50% of each row should be black, and 50% should be white. We also know the bars are only supposed to be 0.5555' pixels wide, so why are we seeing bars that are 5 pixels wide? Well, this is a textbook case of aliasing --- the frequency of 0.9 cycles per pixel is aliased, or folded, around 0.5, which gives us 0.1 (0.5 - (0.9 - 0.5) = 0.1) . We have seen above that 0.1 cycles per pixel gives us 10 pixels per cycle, or bars that are 5 pixels wide.

That explains the width of the bars. But why are they grey, and not black-and-white? The MTF plot of the box filter above provides the answer: at 0.9 cycles per pixel, contrast has dropped to roughly 10%, which manifests as both white and black moving closer to the intensity midpoint of 50% grey (remember, the image shown here are at a gamma value of 2.2). If we see an image such as the one above (box filtered, bar width=0.5555' pixels), how do we know whether we are looking at a pattern with bar width 5 at 10% contrast, or a pattern with bar width 0.5555' pixels at 100% contrast? There is no way of telling, which is why aliasing is a destructive process. By the time you have captured the image with your sensor (without AA filter), you have already lost the ability to tell these two cases apart.

Before concluding, we have to take a quick look at the same pattern, but filtered with a Gaussian point spread function:

Gaussian filtered (standard dev. = 0.568 pixels), bars are 0.5555' pixels wide

If you look closely, you can see that the grey patches are still composed of grey bands of slightly varying shades, but it is immediately clear that the contrast between the bands has dropped significantly. A Gaussian point spread function with a standard deviation of approximately 0.568 has an MTF curve described by the function exp(-6.365x²), which means that at 0.9 cycles per pixel, contrast has dropped to only 0.577%. Incidentally, this Gaussian point spread function appears to be quite close to the AA filter in a Nikon D40 or a Nikon D7000.

Personally, I prefer the fade-to-grey approach. While the box-filtered image (with bars of width 0.5555' pixels) appears to have detail at 0.1 cycles per pixel (bars of width 5 pixels), this detail is false, and was never present in the scene we captured with our sensor.

Conclusion

So which is better: a sensor with an AA filter, or one without? This depends a great deal on what you plan to do with that sensor. A sensor without an AA filter will be more susceptible to aliasing than one with a filter, but this depends critically on the rest of the parameters of the entire optical system.

For example, if your pixel pitch is rather small, say around 1.5 microns, like the sensors found in camera phones, then diffraction will act as a natural AA filter. For a larger pixel pitch of around 5 microns or larger, diffraction will only start acting as an AA filter at small relative apertures (say, beyond f/8, depending on physical sensor size). For large pixels at large apertures, the difference between an AA filter sensor, and one without, boils down to exactly the differences illustrated above. (Update: new post on taking diffraction into account).

Bayer sensors are another matter entirely, but suffice it to say that you must give up some resolution in order to reconstruct colours accurately with this type of sensor.

For photography, I would much rather have a slightly softer image, which can be sharpened afterwards, than a sharper image filled with artifacts that cannot be removed automatically. Remember, blurring the image after capture will not remove the aliasing; the blur must be applied before the sensor samples the incoming scene.

Lastly, always remember that aliasing is only present if your scene contains frequencies above the Nyquist limit of your sensor. Not all images captured without an AA filter will contain aliasing.