Accuracy of sampling images near the Nyquist limit

[alanrockwood]

Oops, my post above go mangled by repeating a quote and repeating part of my response, but hopefully it's not mangled so much as to be impossible to follow.

 

[alanrockwood]

I am posting three images here. The first is if one samples the first function in my first post at the Nyquist rate (a point spacing of 0.125 x_units). The second is if the function is sampled at slightly higher than the Nyquist rate (a point spacing of 0.113636 x_units). The third is the function slightly undersampled (0.1375 x_units).

If you want, you can count the points in each figure to confirm that they are correctly described (sampled at Nyquist limit, sampled slightly above the Nyquist limit, and sampled at slightly below the Nyquist limit.)

As you can see, sampling at slightly above the Nyquist limit gives a beat pattern. So does sampling slightly below the Nyquist limit. Sampling at the Nyquist limit gives an apparently pretty good representation of the function, but that is somewhat deceiving because the phase of the sampling points in this case exactly matched the anti-nodes of the function. If the sampling were shifted slightly it could sample the nodes, and in that case the sampled function is zero at every sampling point, which obviously does not give a good representation of the original function.

The example where the points were sampled slightly above the Nyquist limit would allow one to reconstruct the original function without error if (and only if) a correct reconstruction function is applied.

The example where the points were sampled slightly below the Nyquist limit would not allow one to reconstruct the original function without error because aliasing has occurred. (Actually, exact reconstruction is possible if one has some additional information, but not if the process is done "blind".)

 

[Nodda Duma]

I think the other thing that's bugging me is that your example is not taking into account the effects of an imaging detector on the amplitude of the signal. The fourier transform of a rect function is a sinc function. So, mathematically, you're dropping off the effect that the detector is having on the magnitude of the original pattern. This is a real effect and cannot be ignored like it can (?) for 1-D electronic or audio signals with simple-to-implement downstream filtering.

According to your example, with a Nyquist sampling frequency of 1.1x the original input bandwidth, the maximum amplitude of the sampled signal will only be 9% of the original...if I did my math right (might be off by 2x so peak value of 30% the original maximum amplitude). In a real signal, this contrast is barely detectable by the human eye (threshold of detection is either 10% or 5% depending on who you're talking to). To achieve higher contrast, Nyquist is moved farther out relative to the maximum input bandwidth but your reconstruction of the original frequency gets more accurate, so... You can see this effect in a Moire aka "starburst" pattern.

 

[alanrockwood]

I think the other thing that's bugging me is that your example is not taking into account the effects of an imaging detector on the amplitude of the signal. The fourier transform of a rect function is a sinc function. So, mathematically, you're dropping off the effect that the detector is having on the magnitude of the original pattern. This is a real effect and cannot be ignored like it can (?) for 1-D electronic or audio signals with simple-to-implement downstream filtering.

According to your example, with a Nyquist sampling frequency of 2.2x the original sinusoidal, the maximum amplitude of the sampled signal will only be 9% of the original. In a real signal, this contrast is barely detectable by the human eye (threshold of detection is either 10% or 5% depending on who you're talking to). To achieve higher contrast, Nyquist is moved farther out relative to the maximum input bandwidth but your reconstruction of the original frequency gets more accurate, so... You can see this effect in a Moire aka "starburst" pattern.

No, the Nyquist limit is 2X the frequency of the signal. That is the mathematical definition of the Nyquist limit. The However, the sampling rate was 2.2X the nyquist limit of this signal, i.e. it is sampled at greater than the Nyquist limit, so it satisfies the Nyquist sampling theorem.

 

[alanrockwood]

Regarding the detector response, it is not a sync function. The sync function (sin(x)/x) has both positive and negative values. Neither film nor a digital detector can give a negative response to light. The detector responds as the magnitude (i.e. an absolute value) which has no negative values.

However, the sync function does enter into the reconstruction of the sampled signal. I mentioned an interpolation method that allows the exact reconstruction of the original signal. It uses a brick wall filter, which does involve convolution with a sync function. However, this won't produce negative values in the reconstructed signal, which can only have positive values. One might think that convolution with a sync function would produce some negative values, but they cancel out during the convolution process. (There might be some end-effects in the signal where this might not be strictly true. I haven't investigated the issue of end effects.)

 

[Nodda Duma]

No, the Nyquist limit is 2X the frequency of the signal. That is the mathematical definition of the Nyquist limit. The However, the sampling rate was 2.2X the frequency of the signal, i.e. it is sampled at greater than the Nyquist limit, so it satisfies the Nyquist sampling theorem.

So yeah I was off by 2x. So the output signal magnitude is only 30% of the original. In any case, my point still stands that the real effect is not as accurately reflected in your example.

In addition, you also have to account for the additional contrast-reducing effect of a real optical MTF. This sort of depends on the cutoff frequency of the optics relative to the detector but in practice you see an additional 10-20% reduction in the contrast. So for a real system, the output contrast (because we're talking about an imaging system so lets use the correct terms) is 25%-28% of the magnitude of original sampled scene.

 

[Nodda Duma]

Regarding the detector response, it is not a sync function. The sync function (sin(x)/x) has both positive and negative values. Neither film nor a digital detector can give a negative response to light. The detector responds as the magnitude (i.e. an absolute value) which has no negative values.

Do you really want to go that deep? As I stated in my first post, a forum like this does not do a full discussion justice. Crack open Goodman's book and dive in.

A simple 1-D sinusoidal example is simply not a good representation when talking about imaging systems!

 

[alanrockwood]

So yeah I was off by 2x. So the output signal magnitude is only 30% of the original. In any case, my point still stands that the real effect is not as accurately reflected in your example.

In addition, you also have to account for the additional contrast-reducing effect of a real optical MTF. This sort of depends on the cutoff frequency of the optics relative to the detector but in practice you see an additional 10-20% reduction in the contrast. So for a real system, the output contrast (because we're talking about an imaging system so lets use the correct terms) is 25%-28% of the magnitude of original sampled scene.

I have also not taken into account that a digital imaging detector is not a point sampling device. Each pixel location in the sensor is an area sampling device, which will tend to smooth out some of the effects to some degree, but the effects do not disappear.

 

[alanrockwood]

One thing to be aware of with regard to the beat pattern. The contrast will depend somewhat on how close to the Nyquist limit the sampling is performed. If it is very close to the Nyquist limit there will be large areas in the sampled result were the contrast between sampling points is nearly zero, but in other parts of the beat pattern the contrast between sampling points will be strong.

Of course, as I mentioned in another post, if the image is strongly oversampled then the kind of effects I am bringing up are not very noticable, and one can take the sampled result directly as a quite good representation of the original function, i.e. no fancy reconstruction math is needed to produce a good image result.

As I think you also mentioned in one of your posts, one of the implications of all this is that there is a significant advantage to strongly oversampling the signal. This can reduce aliasing as well as the effect I have been discussing (which is not aliasing, though it might resemble aliasing to some degree.)

Since certain images on film can exceed 100 line pairs per millimeter (by a considerable margin in some cases), a really good scanning system should scan at a rate much greater than 5000 pixels per inch if it is to capture all of the information in the best film images.

I also suspect that the effect I have been discussing here may contribute to the phenomenon sometimes called "grain aliasing". I suspect that grain aliasing may be a result of both true aliasing and the fact that the scanner is trying to record some spatial frequency components that are too close to the Nyquist limit, and hence introduce some distortions into the response. I have not investigated this. It is only a hunch. I have not seen a term for the effect I have brought up in this post. I think a good term might be "pseudo-aliasing" because it mimics aliasing in some respects, but its origin is not from undersampling, so it is not true aliasing.

 

[alanrockwood]

Oh, one other thing. Another implication of this analysis is that it could impact some measurements of resolution that depend on imaging of resolution targets. It would be interesting to analyze this possibility, but I don't want to do it.

 

[wombat2go]

In the interest of collegial discussion let me strongly encourage you to do your own calculation: sample a simple sinusoid at a rate slightly greater than the Nyquist limit, i.e. slightly greater than 2F. When you plot the results I guarantee that you will see a beat pattern.

I had a go at it pictorially. However it took longer than I intended!

Yes I agree that the first function in Alan's #1 post is a RECT() . Cos() ( 1D spatial or time)
To sample it properly would require knowledge of the highest harmonics in that function, not just the Cos().

The sample with green dots in the second plot, is way undersampled, i think, for that function.


I made a 2D equivalent consisting of (4) wavy vertical bars amplitude from peak to zero,
on a pitch of approximately 1/8 of the fundamental .


Then sampling at 1/4000th or a Nyquist of 1/2000th, converted to Fourier.
I can see harmonic content of the function out to about 1300 but presently have no way to measure them.

Here is the reconstituted image with all harmonics
https://app.box.com/s/vm6fqfavk3jr1cmx7ihldan7qkw5fwaf

Then applied a rectangular lo pass filter

This one is passing up to 16th harmonic
https://app.box.com/s/6oclf96ov1m7opg95sjicq4ynlpd0vz4

This one is passing up to 8 harmonics.
https://app.box.com/s/582n8lxy6ckknjtrxg9u2vy7wu8oodxg
8 harmonics is not sufficient to show the wavy bars, but the RECT platform is just there albeit with ramped flanks.

 

[alanrockwood]

I had a go at it pictorially. However it took longer than I intended!

Yes I agree that the first function in Alan's #1 post is a RECT() . Cos() ( 1D spatial or time)
To sample it properly would require knowledge of the highest harmonics in that function, not just the Cos().

The sample with green dots in the second plot, is way undersampled, i think, for that function.


I made a 2D equivalent consisting of (4) wavy vertical bars amplitude from peak to zero,
on a pitch of approximately 1/8 of the fundamental .


Then sampling at 1/4000th or a Nyquist of 1/2000th, converted to Fourier.
I can see harmonic content of the function out to about 1300 but presently have no way to measure them.

Here is the reconstituted image with all harmonics
https://app.box.com/s/vm6fqfavk3jr1cmx7ihldan7qkw5fwaf

Then applied a rectangular lo pass filter

This one is passing up to 16th harmonic
https://app.box.com/s/6oclf96ov1m7opg95sjicq4ynlpd0vz4

This one is passing up to 8 harmonics.
https://app.box.com/s/582n8lxy6ckknjtrxg9u2vy7wu8oodxg
8 harmonics is not sufficient to show the wavy bars, but the RECT platform is just there albeit with ramped flanks.

There are no harmonics in the first graph in my first post. The function plotted in the first figure is Cos(2*pi*F*x)+2, where pi is approximately 3.14159, F is the spatial frequency (4/x_units, just call it 4), and x is the coordinate on the horizontal axis. It is a pure sinusoidal function (offset in the vertical direction by 2 units), so it contains no higher harmonics by definition.

The Nyquist limit for this function is to sample it at a point spacing of 0.125 units in the x direction. (I have been calling the units in the x direction as x_units, but if you want you can just drop the units and call it 0.125.) Any sampling interval smaller than 0.125 satisfies the Nyquist sampling theorem. The second figure is a graph of that function sampled at a point spacing of 0.1136363636. Let's just round that off to 0.1136 for convenience. Since 0.1136 is less than 0.125, the sampling rate in the second graph satisfies the Nyquist theorem.

Go ahead and use a hand calculator to calculate Cos(2*pi*4*x) for x = 0, 0.1136, 0.2272, 0.3408, 0.4544.... Those values for x are the sampling points used for the second figure. (Actually, I used multiples of 0.1136363636, but we are rounding here for convenience.) Then plot the results and you will find that they look just like the second figure in my first post. There is no need to try to work through the problem using graphics. A hand calculator and some graph paper to plot the results will do the job. I used Psi-Plot to do the calculations and to produce the graphs. Excel could probably do it as well. It should only take a few minutes.

 

[wombat2go]

"There are no harmonics in the first graph in my first post. "

It has extensive harmonics becuse it is a rectangular window of 5 periods.

I took liberty to add "zero hold" to your green dots to more closely represent the output from a scanner etc.
https://app.box.com/s/c0vifqhv29m3qhtuagmzfiti85jptjcr

Then ran that into fftw via a front end that takes ASCII files.
Here is the Fourier transform (Amplitude * 500 I think)

https://app.box.com/s/11tydc42l1rsmqacjp8bcnwzvkair4mu

If it is like a continuous radio carrier , then there are finite harmonics.
It was too difficult to extend your 5 cycles to a large number, but I did a similar waveform as a function and verified that
in a window containing 50 cycles, the wanted signal (4 Hz) and the sampler (9 Hz) dominate as diracs.

 

[alanrockwood]

"There are no harmonics in the first graph in my first post. "

It has extensive harmonics becuse it is a rectangular window of 5 periods.

I took liberty to add "zero hold" to your green dots to more closely represent the output from a scanner etc.
https://app.box.com/s/c0vifqhv29m3qhtuagmzfiti85jptjcr

Then ran that into fftw via a front end that takes ASCII files.
Here is the Fourier transform (Amplitude * 500 I think)

https://app.box.com/s/11tydc42l1rsmqacjp8bcnwzvkair4mu

If it is like a continuous radio carrier , then there are finite harmonics.
It was too difficult to extend your 5 cycles to a large number, but I did a similar waveform as a function and verified that
in a window containing 50 cycles, the wanted signal (4 Hz) and the sampler (9 Hz) dominate as diracs.

Wombat,
You are doing the frequency analysis on the wrong figures. None of the figures with the green dots represent the original function. The original function is f(x)=Cos(2*pi*F*x)+2. Five periods of that function are shown in the first first figure.

There are several ways of representing a Fourier series. One is to use Sines and Cosines. All of the frequency components with phase representable by Sine functions in figure 1 have a value of zero. You can take my word for it, or if you want you can calculate the components yourself. This is actually pretty easy to see from symmetry considerations alone, so you don't actually have to apply pencil to paper to calculate those integrals.

The frequency components with phase representable by Cosine functions are given by the following integral.

Where L is the interval over which the analysis is being performed. In this case L=5. If you calculate those integrals you will find that there are only two choices of n for which the results are non-zero. Those are n=0 (which represents the 2 in f(x)) and n=4 (which represents the cos(2*pi*4*x), i.e. the component with a frequency of 4.) Thus, there are neither harmonics nor subharmonics in this function, with the trivial exception of the subharmonic with a frequency of zero.

Now, if you take the sampled results (the ones with the green dots), apply a point-to-point interpolation between the green dots, and then do a frequency analysis on those functions you will get abundant harmonics. You don't even need to calculate the integrals to see that there will be harmonics. All you need to know is that there is an abundance of discontinuities in the slopes, and every discontinuity generates abundant harmonics. But the fact that the sampled results contains abundant harmonics is simply a result of the fact that the sampled results cannot be taken as a direct representation of the original function, which of course is my primary point.

Similarly, if instead of doing a point-to-point interpolation between the dots as your reconstruction, you take each point to represent the center of a "boxcar" element (which is what you have done), then your reconstruction will look like a function with a lot of stairsteps in it. Again, that function has discontinuities that will generate a lot of harmonics, but that is because the reconstruction is not an accurate representation of the original function.

On the other hand, if you do a proper reconstruction, using the sampled function as your input, and provided that you sampled the function above the Nyquist limit, you will re-create the original function, which in this case contains no harmonics. An important point is that you cannot do a proper reconstruction without increasing the number of points in the reconstruction. In fact, rigorously speaking the number of points needs to be infinite, though in practice one can use a smaller number of points to get a quite accurate reconstruction. You can do the reconstruction using the Whittaker-Shannon interpolation formula. This is explained in the Wikipedia article on the Nyquist-Shannon sampling theorem. By the way, if you try to reconstruct the function, but use the original sampling points as your output grid instead of a much larger number of points, you will simply reconstruct the post-sampled function, which is not a good approximation of the original function if the sampling rate is too close to the Nyquist limit.

 

[wombat2go]

Hi Alan,
To:
"You are doing the frequency analysis on the wrong figures. None of the figures with the green dots represent the original function."

In post #36 I did , of a similar but not identical function to your original function,
a 2-D spatial FFT and back again, progressively reducing harmonics, and showing the IFT results.
To:
"if instead of doing a point-to-point interpolation between the dots as your reconstruction,
you take each point to represent the center of a "boxcar" element (which is what you have done),
then your reconstruction will look like a function with a lot of stairsteps in it."

On film "grain", or on a film scanner, or on a cmos camera sensor, present technology, the "no-mans land" between sensor elements is minimised.
The output therefore tends to be more like a classical "zero-hold " sampling with steps.
And not an interpolation between infinitesimal measuring points.

There are artefacts on all technologies, for example "acutance" on b/w film.
Also artefacts on cmos sensor elements of a finite width with small gaps between, as is done presently.

But I don't think any of that relates to your primary point.

Your main drift seems to be:
We can not re build a scanned window,
sampled at 2.2 times the highest wanted frequency in it.
We can re-build an infinitely long train like audio on a radio frequency carrier,
sampled at 2.2 times the highest wanted frequency in it, provided the correct filters are used.

Is that your point?
If so I would agree.

 

[alanrockwood]

Hi Alan,


Your main drift seems to be:
We can not re build a scanned window,
sampled at 2.2 times the highest wanted frequency in it.
We can re-build an infinitely long train like audio on a radio frequency carrier,
sampled at 2.2 times the highest wanted frequency in it, provided the correct filters are used.

Is that your point?
If so I would agree.

No, not quite my main point, though we might be slowly converging toward a meeting of the minds.

You can reconstruct the original function if you sample it at (for example) 2.2 times the highest frequency in the function, but only if you use a proper reconstruction filter. My main point is that it is incorrect to take the sampled result itself and just pretending like it is a correct reconstruction of the original signal. It does not necessarily give a good approximation to the original function, particularly if the sampling rate is only a little bit higher than the Nyquist limit. You have to apply a proper reconstruction filter, and you also have to interpolate points between the sampling points while using the proper reconstruction filter.

This does not invalidate the Nyquist sampling theorem, because all the theorem says is that it is possible to reconstruct the original signal without error, provided that all of the frequency components of the signal are below the Nyquist limit.

So, how does this relate to scanning? Let us suppose that you are photographing a repeating pattern, such as picket fence or a close up of a fabric. Let us suppose the lens is so good that it is not a limiting factor. Let us (for sake of discussion) also assume that there is a perfect brick wall optical filter in front of the sensor, and that the cutoff frequency (spatial frequency) of the filter is at the Nyquist limit for the sensor. If the sensor pitch is just a little better than the Nyquist limit (which of course will depend on how close the camera is to the object.) Your image as it comes off the camera will not produce a good representation of the image because of the effect I have been discussing in this paper. For example, you may find a beat pattern overlaid on the repeated pattern, i.e. a beat pattern that does not exist in the original image. It is pissible to reconstruct the image, but to do so requires applying a proper reconstruction filter and also add interpolating points between the sensor pixels. If you do that (which as far as I am aware is not typically done by image processing software) then you can get a good representation of repetitive features that are close to the frequency of the Nyquist limit (actually, close to half of that frequency.)

In my opinion this calls into question some of the scanner testing that is sometimes reported on the internet. For example, if a set of lines on a resolution target is being imaged by the scanner, and the conditions are close to the Nyquist limit, then beat patterns will appear in the resolution target, and it is only by accident whether one samples in the weak or the strong part of the beat pattern, and this in turn can affect contrast in the image, hence the apparent resolution.

The length of the train is actually not the issue.