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.