Critical angles

It is often said that there is more than one way to skin a cat.

Well, today I discovered an Imatest article that demonstrates just how wildly different slanted edge implementations can (and apparently do) vary. I will leave my critique of said article for another day, but I will note that this article makes reference to the "5 degrees" rule that is often seen when slanted edge measurements are performed.

The "5 degrees" rule states that the orientation of the edge relative to the sensor's photosite grid should be approximately 5 degrees (either horizontal or vertical).

There are two notable reasons for this: firstly, a 5 degree angle is far from the critical angles (the topic of this post), and secondly, a 5 degree angle ensures that the potential non-rotationally symmetric behaviour of the PSF is minimized. A discussion of the non-rotationally symmetric PSFs will also be postponed to a future article.

A closer look at the slanted edge method

Figure 1 illustrates how MTF Mapper constructs the oversampled edge spread function (ESF) that is the starting point of the MTF calculation.

Figure 1: How the ESF is sampled

We want to oversample the ESF so that we can increase the effective Nyquist limit; this is extremely important if we want to measure frequencies close to the natural Nyquist limit of 0.5 cycles per pixel of our sensor. The Shannon-Nyquist theorem shows us that we will have aliasing at frequencies above 0.5 cycles per pixel if we sample at a rate of 1 sample per pixel.

Pushing up our sampling rate to 8x moves the Nyquist limit up to 4 cycles per pixel, which allows us to examine the behaviour of our MTF curve near 0.5 cycles per pixel without fear that we are being misled by aliasing artifacts.

How can we increase the spatial sampling rate of our sensor? Well, we cannot change the sensor, but we can use a trick to generate a synthetic ESF. Looking at Figure 1 above we can see that the edge (represented as a black line) crosses the pixel grid in different places as we move along the edge. More importantly, pay attention to the shortest distance from each black dot (representing the centre of each pixel/photosite) to the black edge. Notice how this distance varies by a fraction of the pixel spacing as we move along the edge.

Let us assume that we have a coordinate system with its origin at the centre of our top/leftmost pixel of our sensor, such that the black dots representing the pixel centres can be addressed by integer coordinates. If we take the (x, y) coordinate of a pixel near the edge, and project this coordinate onto the vector representing the edge normal (i.e., the vector perpendicular to the edge under analysis), then we obtain a real-valued scalar that represents the distance of the pixel centre from our edge. We can pair this projected distance-from-edge value with the intensity of that pixel to form a sample point on our synthetic ESF, as shown in Figure 1.

How does this help us to oversample the ESF? Well, if we choose an appropriate edge orientation angle, say, 5 degrees, then the projected ESF points will be densely spaced. In other words, the average distance between to consecutive samples in our projected ESF will be a fraction of the pixel spacing. We can partition the projected ESF points into bins of width 0.125 pixels to produce a regularly-spaced sampled ESF with 8x oversampling.

We know this works well for 5 degrees (because that is what everyone is doing), but what is so special about 5 degrees? To answer that, we have to slog through some elementary math.

Spacing of projected samples

Figure 2 illustrates one possible way in which we can assign integer coordinates to the pixels near the edge under analysis.

Figure 2: How pixel coordinates are assigned

Note that we can pick an arbitrary origin (shown as (X₀,Y₀) in red); this just simplifies the math that will follow. This point need not fall exactly on the edge, but without loss of generality we can pretend that it does, since this means we can use integer coordinates to refer to the pixel centres of pixels near the edge.

The orientation of the edge can be specified in degrees as measured from the horizontal, but I prefer using the slope of the line. If the angle between the edge and the horizontal is  θ, then the direction perpendicular to the edge can be represented as the unit length vector (-sin(θ), cos(θ)). This would be expressed as a slope 1/Δx = tan(θ), such that Δx = 1/tan(θ).

The normal vector (-sin(θ), cos(θ)) then becomes (-1, Δx) * 1/√(1 + Δx²). We project our pixel centres, represented as integer coordinates (x, y), onto this normal vector by computing the dot product (x, y) · (-1, Δx) * 1/√(1 + Δx²), which evaluates to d(x,y) = 1/√(1 + Δx²) * (-x + yΔx).

The function d(x,y) thus computes the distance that the pixel located at (x, y) is from the origin (X₀,Y₀), which we will pretend falls on the edge; this means that d(x,y) measures the perpendicular distance of point (x, y) from the edge. The projected ESF point is thus [d(x,y), I(x+X₀, y + Y₀)], where I(i, j) denotes the intensity of the pixel located at pixel(i,j).

Suppose that we focus only on the subset of pixels with integer coordinates (p, q) such that 0 ≤ d(p, q) < 1. If we are to achieve 8x oversampling, then there must be at least 8 unique distance values d(p, q) in this interval. In fact, we would require these 8 points to be spread out uniformly such that at least one d(p, q) value falls in the interval [0, 0.125), one in [0.125, 0.25), and so on, such that each of the sub-intervals of length 0.125 between 0 and 1 contain at least one point.

Consider, for example, the case where Δx = 4. This reduces d(p, q) to 1/√(1 + 4²) * (-p + 4q) = (-p + 4q)/√17. Because both p and q are integers, we can deduce that d(p, q) must be an integer multiple of 1/√17. How many integer multiples of 1/√17 can we fit in between 0 and 1? If we enumerate them, we can choose p and q such that (-p + 4q) is the set {0, 1, 2, 3, 4, 5, 6 ...}. But √17 = 4.123106 (and change), so if (-p + 4q) ≥ 5, then d(p, q) > 1. That leaves only the set {0, 1, 2, 3, 4}, such that the only values of 0 ≤ d(p, q) < 1 are {0, 1/√17, 2/√17, 3/√17, 4/√17}.

Whoops! If Δx = 4, then there will only be 5 unique values of d(p, q) between 0 and 1, and we need at least 8 points between 0 and 1 to achieve 8x oversampling! The implications of the failure to achieve 8x oversampling will be covered a bit later; first we must identify the critical angles.

Enumerating the problem angles

We already know that Δx = 4 causes our 8x oversampling to fail; this corresponds to an angle of atan(1/4) = 14.036 degrees. In fact, it is fairly simple to see that for any integer value Δx, we will have Δx + 1 unique values between 0 and 1 (if we include the 0 in our count). For 8x oversampling, the spacing between d(p, q) values must be less than 0.125, which happens when we have at least 8 unique d(p, q) values between 0 and 1. For Δx = 8, we see that 1/√(1 + Δx²) = 1/√65 ≈ 0.12403.

The angles that will lead to a failure of the 8x oversampling mechanism are thus:  45, 26.565051, 18.434949, 14.036243, 11.309932,  9.462322, and  8.130102.

Some other Δx values are also problematic: 1.5, and 2.5. These yield only 2Δx + 1 unique values (including zero). Setting Δx = 1.25 only yields 4Δx + 1 a total of 7 unique values. These fractional slopes occur at angles of 33.69007, 21.80141, and 38.65981 degrees.

There may even be more of these problematic angles, but this is as far as I have come with this analysis. Feel free to comment if you can help me identify other values of Δx that will lead to undersampling.

Dealing with the critical angles

So what exactly happens when we do not have at least one sample every 0.125 pixels along the ESF? The corresponding bin in the resampled ESF will be missing, and leaving gaps in the resampled ESF leads to severe distortion of the MTF because those gaps show up as high-frequency transitions in the FFT.

A workable strategy is to fall back on 4x oversampling. Another strategy is to simply interpolate the from nearby bins. Both of these solutions address the primary issue (gaps in the ESF/PSF), but the residual impact of the interpolation/replacement on the final MTF is harder to mitigate.

A new hope

After my previous post (on improved apodization) I started thinking about the notion of applying low-pass filters to an interpolating function applied directly to the dense ESF samples, before binning is performed. I realized that my explanation of the equivalence between binning and fitting an interpolating function + low-pass filtering + sampling only holds when the points are relatively uniformly distributed within each bin.

This got me thinking that I can probably apply a low-pass filter directly to the dense ESF samples, even before binning. The implementation of this approach feels familiar; it turns out to be similar to the method I implemented to perform importance sampling when using an Airy + photosite aperture PSF (this post). Before describing the new method, first consider this illustration of plain vanilla unweighted binning:

Figure 3: Unweighted binning

The pink boxes denote the bins, each 0.125 pixels wide; the horizontal direction depicted here corresponds to the "d" axis in Figure 2. The midpoint, or representative "x" value for each bin is indicated by the arrows and the values in blue. The green dots represent individual dense ESF samples --- their "y" values are not important in this diagram; the position of the green dots are merely to illustrate where each dense ESF sample is located within each bin in terms of x value, and the number of dots give a rough indication of the density of the dense ESF samples.

If we use plain binning, then we choose as representative x value for each bin the midpoint of the bin. The representative y value is obtained as the mean of the y values of the ESF samples within that bin. In Figure 3, the rightmost bin has many ESF samples quite close to the midpoint of the bin, but almost as many ESF samples near the edge of the bin. The effect of unweighted averaging would be that the samples near the right edge of the bin will carry roughly the same weight as the samples near the middle of our bin, but clearly the samples near the middle of the bin should have had a larger weight in computing the representative value for this bin.

A much better way of binning would be to combine the binning step with the low-pass filtering step. Instead of representing each dense ESF sample as a point, it instead becomes a small rectangle, as shown here:

Figure 4: Weighted binning

Now we can make the weight of each sample point proportional to the overlap of the sample's rectangle and the bin extents. This will allow the samples closer to the midpoint more weight, but it also allows a point to contribute to multiple adjacent bins, depending on the width of the rectangle. This smooths out the transition from one bin to the next, especially if the rectangle is wider than the bin width. (Ok, so the rectangle in the diagram is really just a 1-D interval, not a 2D shape. But the principle still holds.)

Yes, I have just reinvented kernel density estimation. Sigh.

Anyhow, this binning approach also makes the low-pass filtering step explicit, so if each dense ESF sample is now represented by an interval of width w pixels, then we are effectively convolving the ESF with a rect(w * x)  function. We can remove the low-pass filtering effect on the MTF (calculated further down the pipeline) by dividing the MTF by sinc(0.5 * w * f), as I have shown in my previous post.

Our binning process is beginning to look more like a proper approach to sampling: we apply a low-pass filter to our dense ESF points to remove (or at least strongly attenuate) higher frequencies, followed by choosing one representative value at the midpoint of each bin (the downsampling step). By choosing w = 0.33333 pixels, we have a fairly strong low-pass filter, but one that still has a cut-off frequency that is high enough to allow good detail at least up to 3 cycles per pixel.

Because of the (relatively) wide low-pass filter, we could probably drop from 8x oversampling down to 4x oversampling, but I like the extra frequency resolution the 8x oversampling produces in the MTF.

Results

Simulating synthetic images with noise similar to that produced by a D7000 at ISO 800 (but a Gaussian PSF), we can investigate the benefits of the new binning method. Ideally, what we would like to see is no difference between accuracy at a 4 degree angle, and accuracy at one of the critical angles. To quantify this, here is a comparison of 95% percentile of the relative MTF50 error (over a range of MTF50 values from 0.08 cycles/pixel to 0.5 cycles/pixel):

Figure 5: 95% percentile of relative MTF50 error (click to enlarge)

The results are very promising. Most notable is the fact that the new binning method performs virtually identically regardless of edge orientation, with 26.565 degrees being the only angle that is slightly worse than the others.  There may be a slight drop relative to MTF Mapper v0.4.16 (at 4 degrees), but keep in mind the contribution of the change in windowing method discussed in my previous post.

 Just the be sure, I checked for bias at an edge orientation of 4 degrees (although I recycled the ISO800 images):

Figure 6: Relative MTF50 deviation (%)

We can see that the new binning method does not introduce any bias in MTF50 estimates --- of course this is after correction using the MTF of the low-pass filter, as described above.

Conclusion

With the new binning method I can say that MTF Mapper no longer has significant problems with edges of certain orientations. More testing is required, but the 95% percentile of relative MTF50 error appears to be below 5%, regardless of edge orientation, for MTF50 values from 0.08 cycles/pixel through to 0.5 cycles/pixel.

The improved binning method will be included in the next release (which should be v0.4.17).

Improved apodization and bias correction

Following on the relatively recent addition of LSF correction to Imatest, I decided to revisit some of the implementation details of MTF Mapper.

The brutal truth is that MTF Mapper used an empirical correction factor (shock, shock, horror!) to remove the observed bias in measured MTF curves. The empirical correction factor (or rather, family of correction factors) was obtained by generating a synthetic image with a known, analytical MTF curve, and calculating the resulting ratio of the measured curve (as produced by MTF Mapper) to the expected analytical curve.

This had the advantage that it would remove both known distortions, such as that generated by the finite-difference approximation to the derivative (which Imatest refers to as the LSF correction factor), and other distortions which were produced by processes that I did not fully understand at the time.

This post will deal with two of the distortions that I have identified, and I will propose solutions that will enable MTF Mapper to do away with the empirical correction approach.

Apodization

Apodization, also called "windowing", is a way to attenuate some of the artifacts resulting from the application of the FFT (or DFT, if you like) to a signal of a finite length. The DFT/FFT assumes that the signal is periodic, that is, the first (leftmost) sample is preceded (circularly) by the last (rightmost) sample. Applying the FFT to a signal that is discontinuous when treated in this circularly wrapped-around way usually results in significant energy spuriously appearing on the high frequency end of the frequency spectrum.

A common windowing function is the Hamming window, which looks like a cosine function centered on the center of the sequence of samples. The samples are multiplied component-wise with the window function, effectively producing a new set of samples such that the leftmost and rightmost samples are scaled to very low magnitudes. Since the left- and rightmost samples are now all close to zero, we are guaranteed to have a signal that no longer has a discontinuity when wrapping around the left/right ends.

So why would we use apodization as part of the slanted edge method? First, recall how the slanted edge method works:

Step 1: generate the edge spread function (ESF)

This diagram shows how the individual pixel intensities are projected along a line that coincides with the edge we are analyzing. Owing to the angle of the edge relative to the pixel grid, the density of the projected values (along the direction perpendicular to the edge) is much greater than the original pixel spacing. The densely-spaced projected values are binned to form a regularly-spaced set of samples at (usually) 4x or 8x oversampling relative to the pixel grid. This allows us to measure frequencies above the Nyquist limit imposed by the original pixel grid.

Now we can compute the MTF as illustrated here:

Step 2: Compute MTF from PSF using FFT

Notice that the PSF is usually quite compact, i.e., most of the area under the PSF curve is located close to the centre of the PSF curve. This is typical of a PSF extracted from a real-world edge. We see some noise on the tails of the PSF, with visibly more noise on the right side --- this is an artifact of photon shot noise being relative to the signal level, so the noise magnitude is larger in the bright parts of the image.

Anyhow, since the noise is random, we might end up with large values on the edges, such as can be seen on the right end of the PSF samples. This is exactly the scenario which we would like to avoid, so we can apply a window to "squash" the samples near the edges of the PSF.

MTF Mapper had been using a plain Hamming window up to now --- this resulted in a systematic bias in MTF measurements, particularly affecting edges with an MTF50 value below 0.1 cycles per pixel.

Hamming window

Two things are visible here: the noise is suppressed reasonably well (on the ends of the green curve) after multiplying the PSF by the Hamming window function (see right side of illustration), and the PSF appears to contract slightly, effectively becoming slightly narrower after windowing.

The apparent narrowing of the PSF has the expected impact on MTF50 values: they are overestimated slightly.

I identified three possible methods to address this systematic overestimation of MTF50 values (on the low end of MTF50 values): empirical correction (as MTF Mapper has been doing so far), deconvolution, and using a different window function.

We can "reverse" the effect of the windowing after we have applied the FFT to obtain the MTF. By the convolution theorem, we know that convolution in the time domain becomes multiplication in the frequency domain. Since we multiply the PSF by the window function in the time domain, it stands to reason that we must deconvolve the MTF by the Fourier transform of the window function. Except that deconvolution is a black art that is best avoided.

I have tried many different approaches, but the high noise levels in the PSF makes for a poor experience, more apt to inject additional distortion into our MTF than to undo the slight distortion caused by windowing in the first place.


That leaves us only with the last option: choose a different window function. Purely based on aesthetics, I decided on the Tukey window with an alpha parameter of 0.6:

Tukey window

Notice that we may get slightly less noise suppression, but in return we distort the PSF far less. In fact, at this level (MTF50 = 0.05) the distortion is negligible, and no further correction factors are required. This is the new apodization method employed by MTF Mapper.

LSF correction and beyond

As already mentioned, the finite-difference method used to calculate the PSF (or LSF, if you are pedantic) from the ESF is not identical to the ideal analytical derivative of the ESF. A sin(x)/(x) correction factor can be employed to effectively remove this distortion. The Imatest article on this topic does a fine job of explaining the maths behind this correction; the method was originally published by Burns while working at Kodac.

Since MTF Mapper employs 8x oversampling, we must divide the calculated MTF by the function sin(π * f/4)/(π * f/4). Clarification: This stems from the sample spacing that is 0.125 pixels. Plugging this into the finite-difference derivative calculation as explained in the Imatest article we see that for 8x oversampling we will have a correction factor of sin(π * f/4)/(π * f/4) as opposed to the sin(π * f/2)/(π * f/2) we would have had for 4x oversampling.

Even after applying this correction factor, though, we can see a systematic difference between the expected ideal MTF and the MTF produced by the slanted edge method. To understand this (final?) distortion, we have to rewind back to the step where we construct the ESF (helpfully captioned "Step 1" above...).

The projection used to form the dense ESF samples produces a dense set of points, but these points are no longer spaced at convenient regular intervals. The FFT rather depends on being fed regularly spaced samples, so the simplest solution is to bin the samples at our desired oversampling factor. An oversampling factor of 8x thus produces bins that are 0.125 pixels wide.

Again following the path of least resistance, we simply average all the values in each bin to obtain our regularly-sampled ESF. This seems like such a harmless little detail, but if we stop and think about it, we realize that this must be a low-pass filter. Why?

Well, consider first a continuous interpolation function passing through all the ESF samples before binning. We would like to sample this function at regular intervals (0.125 pixels, to be exact), but we know that point sampling will produce horrible aliasing artifacts. The correct approach is to apply a low-pass filter, i.e., convolve our interpolating function with some filter. Let us choose a simple box filter of width 0.125 pixels. If we first convolve the interpolating function with this box filter, and then point-sample at intervals of 0.125 pixels, we end up with exactly the same result as we would obtain from binning followed by averaging all the values in each bin. This approach is optimal in terms of noise suppression for a Gaussian noise source, so even though it sounds simplistic, it is a good solution.

Fortunately, this process is easily reversible by indiscriminate application of the convolution theorem: convolution in the time domain can be reversed by dividing the MTF (in the frequency domain) by the Fourier transform of our low-pass filter. And by now we know that the Fourier transform of a box filter is the sinc() function --- all we have to do is choose the proper frequency.

At 8x oversampling, our bin width is 0.125 pixels, resulting in a low-pass filter of rect(8x). In the Fourier domain, this means we must divide the MTF by sinc(π * f/8) --- this will effectively reverse the attenuation of the MTF induced by the low-pass filter.

To illustrate the effect of these two components (discrete derivative and binning low-pass filter) we can look at a simple example using a Gaussian PSF, with no added noise, and no apoditization. We start with the dense ESF of an edge with an MTF50 value of exactly 0.25 cycles/pixel:

Figure 1: Dense ESF

This ESF is binned into bins of width 0.125 pixels:

Figure 2: binned ESF

Next we calculate the discrete derivative to obtain the PSF:

Figure 3: discrete PSF

 This PSF is passed through the FFT to obtain the following MTF curve:

Figure 4: measured MTF curve

That MTF curve looks pretty good. And it looks very much like half of a Gaussian, just as we would expect. But looks can be deceiving at this scale. We know the true analytical MTF curve that we would expect: a Gaussian with a standard deviation of about 0.2123305 (and change). So next we plot the measured MTF curve divided by the expected MTF curve:

Figure 5: Uncorrected MTF ratio (red)

The dashed blue curve is the sin(π * f/4)/(π * f/4) function, corresponding to the discrete derivative correction, and the red curve is the ratio of measured to expected MTF. Clearly these two curves have roughly the same shape. Let us take our measured MTF curve, divide it by the sinc(f) curve to apply the discrete derivative correction, and plot the ratio of the corrected curve to the expected curve:

Figure 6: Partially corrected MTF ratio (red)

Note how the red curve (corrected MTF divided by expected MTF) has flattened out --- keep in mind that we would expect this curve to flatten out into a straight line. The black dashed line is the function sin(π * f/8)/(π * f/8), i.e., the Fourier transform of the rect(8x) low-pass filter induced by the binning process. Now we can combine the two corrections, i.e., take the measured MTF, divide by the discrete derivative correction, and then divide the result by the low-pass correction; this gives us the "fully corrected" MTF curve. Plotting the fully corrected MTF curve divided by the expected analytical MTF curve yields this:

Figure 7: Fully corrected MTF ratio (red)

The red curve is almost, but not quite, a constant value of 1.0. This demonstrates that the low-pass correction helps to bring us closer to the expected ideal MTF curve.

If we zoom out a bit on the last plot, we see things are not entirely rosy:

Figure 8: Fully corrected MTF ratio (red), wide view

Once we move past a frequency of 1 cycle per pixel, the corrected curve does not match the expected curve so well anymore, at least not when expressed as a ratio. But looking back at Figure 4 above, we see that the measured MTF curve is practically zero beyond 1 cyc/pixel anyway, so we should expect some numerical instability when dividing the measured curve by the expected curve. This also explains my choice of scale in a few of the plots above.

If we express the difference between the fully corrected curve and the expected analytical curve as a percentage of the magnitude of the analytical curve, we see that the fully corrected curve deviates only about 0.15% at 1 cyc/pixel, and only about 0.05% at 0.5 cyc/pixel (Nyquist). For reference, the relative deviation of a completely uncorrected curve is about 10% and 3% at 1 and 0.5 cyc/pixel respectively. Applying only the discrete derivative correction leaves a deviation of about 2.8% and 0.6%.

So adding the correction for the low-pass filter effect of the binning is definitely in the diminishing returns category, but I certainly aim to make MTF Mapper the most accurate tool out there, so no expense is spared.

Summary: The full correction to take care of both the finite-difference correction, and the removal of the attenuation induced by the low-pass filter (implicitly part of the binning operation) is the product of the two individual term, i.e.,

c(f) = sin(π * f/4)/(π * f/4) * sin(π * f/8)/(π * f/8),

The MTF curve is corrected by dividing by this correction factor.

Accuracy evaluation

To demonstrate the effect of the new apoditization and MTF correction approaches, we can look at the MTF50 accuracy over a range of MTF50 values. For each of the MTF50 levels shown below, a number of synthetic images were rendered without adding any simulated noise --- this is to emphasize the inherent bias in measured MTF50 values. All edges were kept at a relative angle of 4.5 degrees, with 30 repetitions rendered using small sub-pixel shifts of the rectangle.

Figure 9: Relative MTF50 deviation on a Gaussian PSF

Our three contestants are MTF Mapper v0.4.16, which employs a Hamming windowing function and empirical MTF curve correction, followed by an implementation that uses a Hamming window with only the discrete derivative correction, and finally the new implementation using a Tukey windowing function with both discrete derivative and binning low-pass corrections.

It is clear that the Hamming window + derivative correction (blue curve) produces a significant bias at low MTF50 values, raising their values artificially (as expected from the apparent narrowing of the PSF). Also note how the MTF50 values are underestimated at higher MTF50 values, which is again consistent with the effects of the binning low-pass filter.

Both the empirical correction method (red curve) and the new Tukey window plus full correction (black curve) display much lower bias in their MTF50 estimates, as seen in Figure 9.

What happens when we use a different PSF to generate our synthetic images? This time I chose the Airy + photosite aperture (square aperture, 100% fill factor) as a representative. This corresponds to something like a D7000 sensor without an OLPF, but without noise.

Figure 10: Relative MTF50 deviation on an Airy+box PSF

Firstly, we see some shockingly large errors on the low MTF50 side. The data points correspond to a simulated aperture of f/64, followed by f/32, f/16, f/8, f/5.6, f/4 and finally f/2.8. A reasonable explanation for the difference between the results in Figure 9 and 10 might be the wider support of the Airy PSF. Typically, the central peak of the Airy PSF is narrower than a Gaussian, but the Gaussian also drops off to zero more quickly, i.e., the Airy PSF has more energy in the tails of the PSF. This means that a wide (f/64) Airy PSF will be affected more strongly by the windowing function, and may even suffer from some truncation of the PSF --- this notion seems to be supported by the difference between the Tukey and Hamming window curves (black vs blue).

Interestingly the empirical correction performed better than expected, doing almost as well as the Tukey + full correction method. This is somewhat unexpected, since the empirical correction factors were calculated from a Gaussian PSF.

Since these experiments were all performed in the absence of simulated noise, they really only test the inherent bias of the various methods. The good news is that the Tukey + full correction approach appears to be an overall improvement over the existing empirical correction, even thought the improvement is really quite small.

Adding in some noise

It always makes sense to look at both bias and variance when comparing the quality of two competing models. In this spirit, the experiments above were repeated under mild noise conditions, corresponding to roughly ISO 800 on a D7000 sensor. First up, the Gaussian PSF:

Figure 11: Standard deviation of relative MTF error on Gaussian PSF

Figure 11 presents the standard deviation of the relative MTF50 error, expressed as a percentage. We see the impact of the Tukey windowing function quite clearly: since the Tukey window does not attenuate such a large part of the PSF (i.e., less of the edge of the PSF is attenuated), we see a small increase in the standard deviation of the relative error. As expected, both the methods using the Hamming window perform nearly identically.

Conclusion

MTF Mapper will employ the new apodization function (Tukey window) as well as the analytically-derived full correction in lieu of the older Hamming window + empirical correction, starting from the next release. This should be v0.4.17 onwards.

The new correction method is more elegant, and makes fewer assumptions regarding the shape of the MTF curve, unlike the empirical correction that was trained on only Gaussian MTFs. But throwing out the empirical correction brings back the strong attenuation of the PSF at lower MTF50 values, so the Hamming window had to be replaced with the Tukey window.

We pay a small price for using the Tukey window, but realistically the MTF50 error should remain below 5% (for an expected MTF50 value of 0.5 c/p) even under quite noisy conditions.

In theory it should be possible to incorporate strong low-pass filtering of the PSF, followed by suitable reversal-via-division of the low-pass filter in the frequency domain. In practice, I have not seen any worthwhile improvement in accuracy. I suspect that some non-linear adaptive filter may be able to strike the right balance, but that will have to wait for now.

Trust, but verify

A while back I wrote:
  "I could not find any synthetic images rendered with specific, exactly known point spread functions. This meant that the only way that I could tell if MTF Mapper was working correctly was to compare its output to other slanted edge implementations." (here)

This was the main motivation for developing mtf_generate_rectangle.

It turns out that Imatest recently updated their SFR measurement algorithm to include the well-known finite-difference-correction (here). I first encountered this correction in one of Burns' papers. According to the Imatest news article, this correction is now included in the ISO 12233:2014 standard.

I have yet to test the new version of Imatest against the synthetic images produced by mtf_generate_rectangle, but I do recall that some quick-and-dirty testing a few years back hinted at Imatest underestimating MTF50 slightly, which would be consistent with an algorithm that does not apply the finite difference correction. As pointed out in the Imatest news article, this difference is really only noticeable when dealing with higher MTF50 values, so this does not imply that all older Imatest results are now suddenly obsolete.


It does raise an important point about traceability and independent verification, though.

Analogue image zoom: comparing 6 megapixel images with 16 megapixel images

The problem

Let us say you have captured a shot of a test chart with a D40, and repeated the process (same lens) with a D7000. The D40 gives you a 6 megapixel image, and the D7000 gives you a 16 megapixel image. You would like to compare the sharpness of the one camera to that of the other.

There several options for executing this comparison:

1. Print both shots at the same size. This is probably the way to go if you intend to print a lot of uncropped photos.
2. Scale down the 16 MP image to 6 MP, and compare at 100% view.
3. Scale up the 6 MP image to 16 MP, and compare at 100% view.
4. Scale both images to some other resolution, e.g., 8 MP (like DxO labs do), or maybe 24 MP.

There is at least one sound reason why option 4 is the better choice amongst options 2 through 4: scaling artifacts. Performing an MTF analysis of an image upscaled with a popular cubic scaling algorithm (Mitchell) reveals that there is some effective contrast enhancement that takes places as part of the scaling process, visible as overshoot and undershoot in an edge profile plot. By scaling both images, you are at least trying to compare apples to apples, especially if you are unsure of exactly what sharpening algorithm your software will employ.

Of course, this entire post deals with visual interpretation of test chart images. If you are interested in other properties (e.g., MTF) then go ahead and use MTF Mapper to perform such measurements directly. The normal slanted edge MTF analysis does not really tell you what your aliasing will look like after dropping the OLPF from your sensor, nor how apparent sharpness is influenced by demosaicing algorithms. For such evaluations visual interpretation might prove useful still.

Another option: simulation

Once you embrace simulated images, such as those produced with mtf_generate_rectangle, you have a fifth option: directly produce a 16 MP (equivalent) image but keep the point spread function (PSF) of the 6 MP camera. I call this feature "analogue scaling", mostly because it effectively resamples the 6 MP image to 16 MP, but without using a discretized 6 MP image as the source. Instead, the usual "infinite precision" analytical description of the target scene is simply scaled down (relative to the photosite pitch), and the sample spacing of the rendered image is adjusted accordingly.

Here is a comparison between a simulated D40 (photosite pitch 7.8 micron, f/4, green light, 4-dot OLPF) and a D7000 (photosite pitch 4.73 micron, f/4, green light, 4-dot OLPF):

Simulated D40 image, effectively magnified by a factor of ~1.65 (click for full-size)

Simulated D7000 image (click for full-size)

These images were generated with the following commands:

mtf_generate_rectangle.exe --target-poly doc/pinch.txt -p airy-box -n 0 --analogue-scale 0.606410 --pixel-pitch 7.8 --aperture 4 -o pinch_d40.png

and

mtf_generate_rectangle.exe --target-poly doc/pinch.txt -p airy-4dot-olpf -n 0  --pixel-pitch 4.73 --aperture 4 -o pinch_d7k.png

respectively. Note that the analogue scale factor is ~0.606, which is 4.73/7.8, i.e., the ratio between the photosite pitch of the D7000 and D40. Specifying an --analogue-scale factor of greater than 1.0 will produce aliasing (and an apparent increase in sharpness), while a factor of less than 1.0 will produce smoothing, as would be expected from upscaling.

Note that the "-n 0" switch turns off simulated sensor noise. Since sensor noise is currently computed in the domain of the output image, this switch is required to produce a correct pair of images. If noise is left on, you would obtain a scaled D40 image with image noise appearing at the size of D7000 pixels, which will clearly cause the D40 image to appear better than it should. This can be fixed (in the mtf_generate_rectangle code) by scaling the "noise image" correctly, but I honestly only thought of this problem now as I am busy writing this blog post :)

Discusssion

Does this approach work? Well, take a look at the point where the converging lines blur into a gray mess in the D40 image (top image above). This appears to happen after the tick marked "2" in the horizontal set of lines --- maybe about one-third of the way from "2" to "1".

In the D7000 image, the extinction point (gray mess) appears almost exactly at the tick marked "1" in the horizontal set of converging lines. This appears about right, since we know the linear resolution of the D40 is about 0.606 that of the D7000. Not really a rigorous proof, but at least reassuring.

To summarize: the --analogue-scale option of mtf_generate_rectangle will allow you to produce a pair of synthetic images of a test chart to simulate two different cameras, with different photosite pitch values, but without the hassle or potential artifacts introduced by upscaling the discrete image produced by the lower-resolution camera.

Of course, this type of simulation will allow you to investigate potential future sensors too. How would a 50 MP APS-C camera render a resolution test chart .... ?

ps: MTF Mapper version 0.4.16 is finally available here --- this is the first version to support all the required features to produce the results found in this post.

mtf_generate_rectangle grows up

Fed up with squares?

If you have used mtf_generate_rectangle before you will know why it is called mtf_generate_rectangle. A bit over two years ago, when I started working on the first release of MTF Mapper, I ran into a specific problem: I could not find any synthetic images rendered with specific, exactly known point spread functions. This meant that the only way that I could tell if MTF Mapper was working correctly was to compare its output to other slanted edge implementations.

While this is sufficient for some, it did not sit well with me. What if all those other implementations were tested in the same way? If that was the case, then all slanted edge implementations (available on the Internet) could be flawed. Clearly, some means of verifying MTF Mapper independently of other slanted edge algorithm implementations was required.

Thus mtf_generate_rectangle was born. The original implementation relied on sampling a dense grid, but quickly progressed to an importance sampling rendering algorithm tailored to a Gaussian PSF. The Gaussian PSF had some important advantages over other (perhaps more realistic) PSFs: the analytical MTF curve was know, and simple to compute (the MTF of a Gaussian PSF is just another scaled Gaussian). Since the slanted edge algorithm requires a step edge as input, it seemed logical to choose a square as the target shape; this would give us four step edges for the price of one.

Such synthetic images, composed of a black square on a white background, are perfectly suited to the task of testing the slanted edge algorithm. Unfortunately, they are not great for illustrating the visual differences between different PSFs. There are a few well-know target patterns that are found on resolution test charts designed for visual interpretation. The USAF1951 pattern consists of sets of three short bars (see examples later on); the width of these bars are decreased in a geometric progression, and the user is supposed to note the scale at which the bars are no longer clearly distinguishable.

Another popular test pattern is the Siemens star. This pattern comprises circular wedges radiating from the centre of the design. The main advantage of the Siemens star is that resolution (spacing between the edges of the wedges) decreases in a continuous fashion, as opposed to the discrete intervals of the USAF1951 chart. I am not a huge fan of the Siemens star, though, mostly because it is hard to determine the exact point at which the converging bars (wedges) blur into a gray mess. It is far too easy to confuse aliasing with real detail on this type of chart. Nevertheless, other people seem to like this chart.

Lastly, there is the familiar "pinched wedge" pattern (also illustrated later in this post), which contains a set of asymptotically convergent bars. The rate of convergence is much slower than the Siemens star, and a resolution scale usually accompanies the pattern, making it possible to visually measure resolution in a fashion similar to the USAF1951 chart, but with slightly more accuracy. I rather like this design, if only for the fact that the resulting pictures are aesthetically pleasing.

Today I announce the introduction of a fairly powerful new feature in mtf_generate_rectangle: the ability to render arbitrary polygon shapes.

The implementation

You can now pass fairly general polygonal targets using the "--target_poly <filename>" command line option. The format of the file specified using <filename> is straightforward: any number of polygons can be specified, with each polygon defined by an integer <n> denoting the number of vertices, followed by <n> pairs of <x,y> coordinates. I usually separate these components with newline characters, but this is not critical.

The polygons themselves can be convex or concave. In theory, non-trivial self intersections should be supported, but I have not tested this myself yet. There is currently no way to associate multiple contours with a single polygon, thus you cannot specify any polygons containing holes. I work around this by simply splitting the input polygons with a line passing through such a hole: for example, a polygon representing the number "0" can be split down the middle, producing two concave polygons that touch on the split line.

General polygon intersections

For a while now I have cheated by relying on the Sutherland-Hodgeman algorithm to compute the intersection between two polygons. Specifically, this operation is required by all the importance sampling algorithms involving a non-degenerate photosite aperture (e.g., "airy-box" and "airy-4dot-olpf" PSF options specified with the "-p" command to mtf_generate_rectangle). This article explains the process in more detail, but the gist is that each "sample" during the rendering process is proportional to the area of the intersection between the target polygon geometry and a polygon representing the photosite aperture (suitably translated). If we assume that the photosite aperture polygon is simply a square (or more general, convex) then we can rely on the Sutherland-Hodgeman algorithm to compute the intersection: we simply "clip" the target polygon with the photosite aperture polygon, and compute the area of the resulting clipped polygon.

Now this is where the cheat comes in: the clipped result produced by the Sutherland-Hodgeman algorithm is only correct if both polygons are convex. If the target polygon (the clipee) is concave, and the clipping polygon is convex, the Sutherland-Hodgeman algorithm may produce degenerate vertices. (see figure 1 below). The cheat that I employed relied on the observation that degenerate sections of a polygon have zero area, thus they have no influence on the sampling process.

Fig 1: Clipping a concave polygon with the Sutherland-Hodgeman algorithm

This allowed me to continue using the Sutherland-Hodgeman algorithm, mostly because the algorithm is simple to implement, and very efficient. This efficiency stems from the way in which an actual intersection vertex is only computed if an edge of the clipee polygon is known to cross an edge of the clipper polygon. This is a huge saving, especially if one considers that the photosite aperture polygon (clipper) will typically be either entirely outside the target polygon, or entirely inside the target polygon, except of course near the edges of the target polygon.

All this comes apart when the photosite aperture polygon becomes concave. To solve that problem requires significantly more effort. For a start, a more general polygon clipping algorithm is required. I chose the Greiner-Hormann algorithm, with Kim's extensions to cater for the degenerate cases. In this context, the degenerate cases occur when some of the vertices of the clipee polygon coincide with vertices (or edges) of the clipper polygon. This happens fairly often when constructing a quadtree (more on that later). At any rate, the original Greiner-Hormann algorithm is fairly straightforward to implement, but adding Kim's enhancements for handling the degenerate cases required a substantial amount of effort (read: hours of debugging). The Greiner-Hormann algorithm is quite elegant, and I can highly recommend reading the original paper.

Internally, mtf_generate_rectangle classifies polygons as being either convex or concave. If the photosite aperture is convex, the Sutherland-Hodgeman algorithm is employed during the sampling process, otherwise it will fall back to Kim's version of the Greiner-Hormann algorithm. The performance impact is significant: concave photosite polygons render 20 times more slowly than square photosite polygons when rendering complex scenes. For simpler scenes, the concave photosite polygons will render about four times more slowly than squares; circular (well, 60-sided regular polygons, actually) will render about two times more slowly than squares.

Part of this difference is due to the asymptotic complexity of the two clipping algorithms, expressed in terms of the number of intersection point calculations: the Sutherland-Hodgeman algorithm has a complexity of O(c*n), where "c" is the number of crossing edges, i.e., c << m, where "m" is the number of vertices in the clipee polygon, and "n" is the number of edges in the clipper. The Greiner-Hormann algorithm has a complexity of O(n*m); on top of that, each intersection vertex requires a significant amount of additional processing.

Divide and conquer

To offset some of the additional complexity of allowing arbitrary target polygons to be specified, a quadtree spatial index was introduced. The quadtree does for 2D searches what a binary tree does for linear searches: it reduces the number of operations from O(n) to O(log(n)).

First up, each polygon is wrapped with an axis-aligned bounding box (AABB), which is just an educated-sounding way of saying that the minimum and maximum values of the vertices are recorded for both x and y dimensions of a polygon. This step already offers us a tremendous potential speed-up, because two polygons can only intersect if their bounds overlap. The bounds check is reduced to four comparisons, which can be implemented using short-circuit boolean operators, so non-overlapping bounds can be detected with as little as a single comparison in the best case.

Once each individual polygon has a bounding box, we can start to aggregate them into a scene (internally, mtf_generate_rectangle treats this as a multipolygon with its own bounding box). The quadtree algorithm starts with this global bounding box, and splits it into four quadrants. The bounding box of each quadrant is taken as a clipping polygon, clipping all the polygons to fit exactly inside the quadrant.

After one iteration, we have potentially reduced the number of intersection tests by a factor of four. For example, if we determine (using the bounding boxes) that the photosite aperture polygon falls entirely inside the top-right quadrant, then we only have to process the (clipped) polygons found inside that quadrant. If a quadrant is empty, we can simply skip it; otherwise, we can shrink the bounding box to fit tightly around the remaining clipped polygons (see figure 2 below).

The next logical step is to apply this quadrant subdivision recursively to each of the original quadrants. We can keep on recursively subdividing the quadrants until a manageable number of polygons (or more correctly, a manageable number of polygon edges) is reached in each quadrant subtree. We must balance the cost of further subdivision against the gains of reducing the number of edges in each subdivided quadrant. Every time that we add another level to the quadtree we add four additional bounds checks --- eventually the cost of the bounds checks add up.

Fig 2: Levels 0 (green), 1 (blue) and 2 (magenta) of the Quadtree decomposition of the scene (light gray)

Fig 3: Levels 2 (green), 5 (blue), and 6 (magenta) of the Quadtree decomposition

If the number of quadtree levels is sufficient, we end up "covering" the target polygons with rectangular tiles (the bounding boxes of the quadrants), providing a coarse approximation of the target polygon shape (see figure 3 above). Every sampling location outside these bounding boxes can be discarded early on, so rendering time is not really influenced by the size of the "background" any more.

If the quadtree is well-balanced, the amount of work (number of actual polygon-polygon intersection tests) can be kept almost constant throughout the entire rendered image, regardless of the number of vertices in the scene. I have confirmed this with some quick-and-dirty tests: halving the number of vertices in a scene (by using a polygon simplification method) has almost no impact on rendering time.

Some examples

Enough talk. Time for some images:

USAF1951-style chart with 4 levels

Siemens star chart

Pinch chart (click for full size)

All the charts above were rendered with a photosite pitch of 4.73 micron, using the Airy + square photosite (100% fill-factor) model at an aperture of f/4, simulating light at 550 nm wavelength. The command for generating the last chart would look something like this:

mtf_generate_rectangle.exe --pixel-pitch 4.73 --aperture 4 -p airy-box -n 0 --target-poly pinch.txt 

where "pinch.txt" is the file specifying the polygon geometry (which happens to be in the same folder as mtf_generate_rectangle.exe in my example above). These target polygon geometry files are included in the MTF Mapper package from version 0.4.16 onwards (their names are usaf1951r.txt, siemens.txt, and pinch.txt.)

OLPF demonstration

The "pinch chart" provides a very clear demonstration of the effects of the Optical Low-Pass Filter (OLPF) found on many DSLRs (actually, most, depending on when you read this).

Before I present some images, first a note about effective chart magnification. Most real-world test charts are printed at a known size, i.e., you can say with confidence that a particular target (say, a USAF1951 pattern) has a known physical size, and thus measures physical resolution expressed in line pairs per millimetre (lp/mm). It is relatively straightforward to extend this to synthetic images generated with mtf_generate_rectangle by carefully scaling your target polygon dimensions. For the time being, though, I prefer to fall back to a pixel-centric view of the universe. In other words, I chose to specify my target polygon geometry in terms of pixel dimensions. This was mostly motivated by my desire to illustrate specific effects (aliasing, etc.) visually. Just keep that in mind: the images I present below are not intended to express resolution in physical units; they are pretty pictures.

With that out of the way, here is a sample of a hypothetical D7000 without an OLPF --- this could be something like the Pentax k5-IIs.

Hypothetical D7000 without OLPF, f/4

And here is the same simulated image, but this time using an OLPF, i.e., this should be quite close to the real D7000:

D7000 (with OLPF), f/4

I repeated the simulations using an f/8 aperture:

Hypothetical D7000 without OLPF, f/8

and again, the D7000 (with OLPF) simulated at f/8:

D7000 (with OLPF), f/8

Here is a crop comparing the interesting part of the chart across these four configurations:

First up, notice the false detail in the "f/4, no OLPF" panel, occurring to the right of the scale bar tick labelled "1". This is a good example of aliasing --- compare that to the "f/4, OLPF" panel, which just fades to gray mush to the right of its tick mark. In the bottom two panels we can see the situation is significantly improved at f/8, where diffraction suppresses most of the objectionable aliasing.

Simulating microlenses: kicking it up a notch

Preamble

My first stab at simulating microlenses made some strong assumptions regarding the effective shape of the photosite aperture. Reader IlliasG subsequently pointed me to an illustration depicting a more realistic photosite aperture shape --- which happens to be a concave polygon.

At first, it might seem trivial to use this photosite aperture shape in the usual importance sampling algorithm employed by mtf_generate_rectangle. It turns out to be a bit more involved than that ....

The mtf_generate_rectangle tool relied on an implementation of the Sutherland-Hodgeman polygon clipping routine to compute the area of the intersection of the photosite aperture and the target polygon (which is typically a rectangle). The Sutherland-Hodgeman algorithm is simple to implement, and reasonably efficient, but it requires the clipping polygon to be convex, so I required a new polygon clipping routine to allow concave/concave polygon intersections (astute readers may spot that I could simply exchange the clipping/clippee polygons, but I wanted concave/concave intersections anyway). After some reading, it seemed that the Greiner-Hormann algorithm had a fairly simple implementation ...

... but it did not handle the degenerate cases (vertices of clipping/clippee polygons coinciding, or a vertex falling on the edge of the other polygon). Kim's extension solves that problem, but it took me a while to implement.

Effective photosite aperture (with microlenses)

The Suede (on dpreview forums) posted a diagram of the effective aperture shape after taking the microlenses into account. I thumb-sucked an analytical form for this shape, which looks like this (my shape in cyan overlaid on The Suede's image):

The fit of my thumb-sucked approximation is not perfect, but I declare it to be good enough for government work. I decided to call this the rounded-square photosite aperture (that is the identifier used by mtf_generate_rectangle).

I am not sure how to scale this shape relative to the 100% fill-factor square. Intuitively, it seems that the shape should remain inscribed within the square photosite, or otherwise the microlens would be collecting light from the neighbouring photosites too. This type of scaling (as illustrated above) still leaves the corners of the photosite somewhat darkened, which is what we were aiming for. Incidentally, this scaling only gives me a fill-factor of ~89.5%. I guess the "100% fill-factor" claim sometimes seen in connection with microlenses applies to equivalent light-gathering ability, rather than geometric area.

Results

MTF curves for 0-degree step edge

MTF curves for 45-degree step edge

The two plots above illustrate the MTF curves of three possible photosite aperture shapes, combined with an Airy PSF (aperture=f/5.6, photosite pitch=4.73 micron, lambda=550 nm). The first plot is obtained by orienting the step edge at 0 degrees, i.e., our MTF cross-section is along the x-axis of the photosite. In the second plot, the step edge was oriented at 45 degrees relative the photosite, i.e., it represents the diagonal across the photosite.
Both plots include the MTF curves for an inscribed circle aperture, for comparison. Note that the fill-factors have not been normalized, that is, each aperture appears at its native size, which maximizes aperture area without going outside the square photosite's bounds.

Purely based on its fill factor of ~90%, we would expect the first zero of the rounded-square aperture's MTF curve to land between the 100% fill-factor square and the 78% fill factor circle, which is clearly visible in the first plot. In fact, the rounded-square aperture's MTF curve appears to be a blend of the square and circle curves, which makes sense.

The second plot above shows that the rounded-square aperture still exhibits some anisotropic behaviour, but that the effect is less pronounced than that observed with a square photosite (see this article for more details on anisotropic behaviour); this also seems logical given the shape.

In the real world (well, simulated real world, at least)

The MTF curves show some small but measurable differences between the 100% fill-factor square photosite aperture and the ~90% rounded-square photosite aperture response to a step edge. But can you see these differences in an image?

USAF1951-style chart, f/1.4, 100% fill-factor square photosite aperture

USAF1951-style chart, f/1.4, ~90% fill-factor rounded square photosite aperture

Well ... not really (click on the image to see full-size version). I even opened up the aperture to f/1.4 to accentuate the differences in the photosite apertures. Just to show you something, here is a rendering using a highly astigmatic photosite aperture (a rectangle that is 0.01 times the photosite pitch in height, but one times the pitch wide):

USAF1951-style chart, f/1.4, 2% fill-factor thin rectangular photosite aperture

Note that this is basically point-sampling in the vertical direction, but box-sampling in the horizontal direction. This shows up as rather severe aliasing (jaggies) in the vertical direction.

In the real real world

So how do these simulated MTF curves compare to actual measured MTF curves? In a previous post I described the method I used to capture the MTF of a Nikon D40 camera with a sharp lens set to an aperture of f/4. Here is a comparison of the simulated MTF curves to the empirically measured MTF curve.

Nikon D40 at f/4

At first glance it might seem that the 100% fill-factor square photosite aperture simulation is marginally closer to the measured curve, but keep in mind that these simulations were both performed with an OLPF split factor of 0.375. This value of 0.375 was determined by trial and error using the 100% fill-factor photosite simulation --- it is likely that the optimal OLPF split factor for the rounded-square photosite aperture model is different. In fact, I would expect a slightly larger value, say around 0.38 or 0.385 to perform better, purely on the difference in fill factor (100% vs ~90%) between the two simulations.

So yes, you could say I am lazy for not optimizing the OLPF split factor for the rounded-square photosite aperture model right now, but I do not feel comfortable doing any sort of quantitative comparison between the models with only one empirical sample at hand (one measured D40 MTF curve). Until such time as I have sufficient data to perform a proper optimization and evaluation of the models, I will leave it at the following statement: it certainly appears that the rounded-square model is a viable approximation of the photosite aperture of the D40.