Super-resolution imaging using the spatial-frequency filtered intensity fluctuation correlation.

We report an experimental demonstration of a nonclassical imaging mechanism with super-resolving power beyond the Rayleigh limit. When the classical image is completely blurred out due to the use of a small imaging lens, by taking advantage of the intensity fluctuation correlation of thermal light, the demonstrated camera recovered the image of the resolution testing gauge. This method could be adapted to long distance imaging, such as satellite imaging, which requires large diameter camera lenses to achieve high image resolution.

Improving the resolution of optical imaging has been a popular research topic in recent years123456. A commonly used simple approach is to measure the autocorrelation of two identical classical images, effectively squaring the classical image, 〈I1(1)〉 〈I1(1)〉, where 1 is the transverse coordinate of the detector. This autocorrelation produces a maximum gain of the spatial resolution. However, the imaging resolution of such a setup can be further improved by changing the measurement from 〈I1(1)〉 〈I1(1)〉, in terms of intensity, or 〈n1(1)〉 〈n1(1)〉, in terms of photon number counting, to the intensity fluctuation correlation 〈ΔI1(1)ΔI2(2)〉, or 〈Δn1(1)Δn2(2)〉, where 1 and 2 are the transverse coordinates of two spatially separated detectors. Then, if only those fluctuation correlations due to the higher spatial frequencies from ΔI2(2) are selected, a super-resolving image can be observed from the joint detection of the intensity fluctuations at the two detectors. The physics behind this super-resolution is similar to the original thermal light ghost imaging78, and is quite different from an autocorrelation measurement. It should be emphasized that the reported result is also different than that of Oh et al.6; while the authors measured the intensity fluctuations, it was the intensity fluctuation autocorrelation 〈ΔI1(1)2〉, which was still limited by the resolution improvement of an autocorrelation measurement.

In this Report, we demonstrate a camera with resolution beyond the classical Rayleigh limit. Similar to the original thermal ghost imaging experiments78, the camera produces an image by the measurement of 〈ΔI1(1)ΔI2(2)〉; the camera consists of a typical imaging setup, except it has two sets of independent and spatially separated detectors: D1 placed on the image plane, and D2 placed on the Fourier transform plane. Crucially, D2 integrates (sums) only the higher spatial frequencies, or transverse wavevectors, by blocking the central area of the Fourier transform plane. The image is calculated from the intensity fluctuations of D1, at each transverse position 1, and the bucket detector D2. The measurement can be formulated as ΔR(1) = 〈ΔI1(1) ∫d2F(2)ΔI2(2)〉, where F(2) is a filter function which selects the higher spatial frequencies.

Experimental setup

Figure 1 illustrates the laboratory-demonstrated camera setup. We used a standard narrow spectral bandwidth pseudo-thermal light source consisting of a 10 mm diameter 532 nm wavelength laser beam scattered by millions of tiny diffusers on the surface of a rotating ground glass. The object imaged was the 5–3 element of a 1951 USAF Resolution Testing Gauge. Like a traditional camera, the imaging lens, L, had an aperture limited by an adjustable pinhole to approximately 1.36 mm diameter. However, this imaging device has two optical arms behind its imaging lens L. The light transmitted by the object falls on a single-mode 0.005 mm diameter fiber tip that scans the image plane of arm one and is then interfaced with a photon counting detector D1. The combination of the scanning fiber tip, D1, and photon counting detector acts as a CCD array. In arm two, the second lens L is placed behind the image plane and performs a Fourier transform of the field distribution of the image plane of arm two. D2, a scannable multimode 0.105 mm diameter fiber interfaced with a photon counting detector, is placed in the Fourier transform plane of L and integrates only the higher spatial frequencies while filtering out the lower spatial frequency modes; we emphasize that the placement of this “spatial filter” does not depend on knowledge of the Fourier transform of the image.

Figure 1

Experimental setup: a 10 mm diameter 532 nm wavelength laser beam scatters from a rotating ground glass (i) and strikes a 1951 USAF Resolution Testing Gauge (ii), then the imaging lens L (iii) and a pinhole of ≈1.36 mm diameter (iv).The light splits to arm one, where it is collected by a scanning point detector D1, and arm two, where it is collected by a spatially filtered bucket detector, consisting of a lens L located on the image plane and a multimode fiber tip in the Fourier transform plane. The image is then calculated using the Photon Number Fluctuation Correlation (PNFC) circuit.

The intensity fluctuations in this experiment were recorded by a Photon Number Fluctuation Correlation (PNFC) circuit910, which independently records the arrival time of each photo-detection event at D1 or D2. The intensity, measured by the number of photons detected per second, is divided into a sequence of short time windows, each of which needs to be less than the second-order coherence time of the light; it is important that the width of the time window not be too long. The software first calculates the average intensity per short time window, , where s = 1, 2 indicates the detector, and then the difference or fluctuation term for each time window: . The corresponding statistical average of 〈Δn1Δn2〉 is thus

It should be emphasized that when the fluctuation correlation is calculated between the two detectors, Δn1 and Δn2 have not yet been time-averaged. The time averaging is performed after the correlation, as indicated by the multiplication appearing inside the sum over j.

Experimental results

Typical experimental results are presented in Fig. 2. In this measurement, the 5–3 element of a 1951 USAF Resolution Test Target was imaged in one dimension by scanning D1 in the x-direction along the slits. Figure 2(a) shows a completely unresolved classical image of the three slits, I1(x1), that was directly measured by the scanning detector D1. For reference the gray shading indicates the location of the slits. Figure 2(b) shows two results: the black dots plot the autocorrelation [I1(x1)]2, while the blue triangles show the fluctuation autocorrelation 〈[ΔI1(x1)]2〉 at each point. 〈[ΔI1(x1)]2〉 was calculated from the intensity fluctuation autocorrelation of D1. The measurements in Fig. 2(b) have a resolution gain, similar to that of Oh et al.6. Using the Rayleigh limit1112 δx = 1.22λs/D1112, where s is the distance from lens L to the image plane and λ is the wavelength of illumination, the expected resolution of the autocorrelation in the image plane is approximately mm which, as seen in 2(b), is not enough to resolve the three slits which have a slit-to-slit separation of about 0.13 mm. However, by spatially filtering arm two, the three slits of the 5–3 element of the gauge were clearly separated when correlated with arm 1, as seen in Fig. 2(c). The error bars in Fig. 2(a) and (b) are quite small, especially when compared to those in Fig. 2(c). This is a typical negative feature of second-order measurements; compared to first-order classical imaging, in order to achieve the same level of statistics the reported imaging mechanism needs a longer exposure time. How much longer the measurement takes depends on the power of the light source and other experimental parameters.

Figure 2

Resolution comparison for different imaging methods of three 0.01241 mm wide slits imaged by a 10 mm diameter source: (a) unresolved first-order classical image, where the gray shading marks the location of the slits; (b) unresolved images from the fluctuation autocorrelation. The black dots indicates [I1(x1)]2 and the blue triangles show [ΔI1(x1)]2, as seen in Oh et al.6; (c) completely resolved image observed from ΔR(x1) where the solid line is a Gaussian data fitting.

Figure 2(c) is the sum of two measurements obtained by placing the bucket fiber tip at two points in the Fourier transform plane: x2+ = 0.22 mm and x2− = −0.24 mm; this selects the higher spatial frequency modes which fall onto the two fiber tips and “blocks” all other spatial frequency modes. We represent this mathematically with the filter function F(x2) = Π(x2 − x2+, D) + Π(x2 − x2−, D) to simulate the physical “spatial filtering”, where Π(x, w) is a rectangle function of width w, D = 0.105 mm is the fiber diameter, and x2 is measured from the central maximum. Then in one dimension for  ∝ kx2/f, where f is the focal length of the bucket lens, ΔR(x1) = 〈ΔI1(x1)∫dx2F(x2)ΔI2(x2)〉, and only integrates the higher spatial frequencies collected by the bucket detector in the neighborhood of x2 = x2+, x2−. Again, this “spatial filtering” does not require any knowledge of the object or its Fourier transform function.

One way to improve these results is to replace D2 with a CCD array; the CCD would still be in the Fourier transform plane, but with the central pixels blocked. This would allow D2 to collect more light of higher spatial frequencies. Although the limits of our equipment, software data storage, and time constraints prevented the authors from making such improvements, in this reported measurement all three slits of the resolution gauge are certainly well-resolved, while both the classical imaging and the autocorrelation mechanisms could not resolve it.

Discussion and theory

In the experiment, we use the spatial correlation of the noise, 〈ΔI1(r1, t1)ΔI2 (r2, t2)〉, to produce an image from the joint photo-detection of two independent and spatially separated photodetectors, D1 and D2. In the following, we outline the theory behind our experiment. First we briefly consider how a first-order or classical camera produces an image in its image plane, 〈I1(1)〉.

The experiment was performed using a pseudothermal light source, created by placing a rotating ground glass in the path of a laser beam. The ground glass contains a large number of tiny scattering diffusers which act as sub-sources. The wavepackets of scattered light play the role of subfields; each diffuser scatters a subfield to all possible directions, resulting in the subfields acquiring random phases. Each sub-field propagates from the source plane to the image plane by means of a propagator or Green’s function, E(1) = Eg(1), where E is the initial phase and amplitude of the field emitted by sub-source m and g(1) is the Green’s function which propagates the light from the m sub-source located at to the point 1 at some distance z from the source plane. To simplify the problem, we assume the fields are monochromatic and ignore the temporal part of the propagator. Then the light measured at coordinate (r, t) is the result of the superposition of a large number of subfields, , each emitted from a point sub-source,

where 〈I(r, t)〉, the mean intensity, is the result of the mth subfield interfering with itself; ΔI(r, t), the intensity fluctuation, is the result of the mth subfield interfering with the nth subfield, m ≠ n, and is usually considered noise because 〈ΔI(r, t)〉 = 0 when taking into account all possible random phases of the subfields.

A classical imaging system measures the mean intensity distribution on the image plane, 〈I1(1)〉, where we have assumed a point detector D1 is placed at coordinate 1, the transverse coordinate of the image plane. In an ideal imaging system, the self-interference of subfields produces a perfect point-to-point image-forming function. The ideal classical image assuming an infinite lens is the convolution between the aperture function of the object |A()|2 and the image-forming δ-function which characterizes the point-to-point relationship between the object plane and the image plane111213.

where μ = s/s is the magnification factor, g() is a Green’s function propagating the m subfield from the source plane to the object plane over a distance z. and g(1) is a function propagating the subfield from the object plane to the detection plane over a distance s + s, and including the imaging lens. A() is an arbitrary function describing the object aperture.

In reality, due to the finite size of the imaging system, we rarely have a perfect point-to-point correspondence. Incomplete constructive-destructive interference blurs the point-to-point correspondence to point-to-spot correspondence. The δ-function in the convolution of Eq. 3 is then replaced by a point-to-spot image-forming function, or a point-spread function which is determined by the shape and size of the lens. For a lens with a finite diameter, one common model describes the shape or pupil of the lens as a disk of diameter D:

where the sombrero-like point-spread function is defined as somb(x) ≡ 2J1(x)/x; J1(x) is the first-order Bessel function. The image resolution is determined by the width of the somb-function: the narrower the higher. A larger diameter lens results in a narrower somb-function and thus produces images with higher spatial resolution.

To simplify the mathematics, it is common to approximate a finite lens as a Gaussian with diameter D, but a smoother falloff than the disk approximation. This leads to a Gaussian imaging-forming function:

This Gaussian version of the imaging equation will be used later in numerical calculations to simplify the mathematical evaluation.

It is clear from Eqs 4 and 5 that for a chosen value of distance s, a larger imaging lens and shorter wavelength will result in a narrower point-spread function, and thus a higher spatial resolution of the image.

Now we consider the noise produced image that is observed from the measurement of Fig. 1 by means of 〈ΔI1(1)ΔI2(2)〉. To make the explanation of the experimental results easier to follow, first we examine the case where two point scanning detectors D1 and D2 are placed in the image planes of arm one and arm two:

The calculation of is straightforward:

Next, we complete the summation over m in terms of the subfields, or the sub-sources, by means of an integral over the entire source plane. This integral results in the well-known Hanbury-Brown Twiss (HBT) correlation: somb2[(πΔθ)/λ| − |], where Δθ is the angular diameter of the light source relative to the object plane. To simplify further calculations, we assume a large value of Δθ and approximate the somb-function to a δ-function evaluated at  = ,  = ′. 〈ΔI1(1)ΔI2(2)〉 is therefore approximately equal to:

It is clear that when 1 = 2 in Eq. 8, the measurement of 〈ΔI1(1)ΔI2(1)〉 produces an image with a resolution gain, with an imaging resolution due to the image-forming somb-functions, i.e., . When the lens is large enough to resolve the object, the result is a point-to-point reproduction of the image only when 1 = 2; otherwise for small lens apertures Eq. 8 forms a point-to-spot image when | + /μ| < λs/D and | + /μ| < λs/D.

Now we move D2 to the Fourier transform plane of L of arm two, i.e., to its focal plane, effectively performing a Fourier transform of the field distribution of the image plane. In addition, D2 is placed off-center relative to the optic axis of the lens to select part of the spatial frequencies on the Fourier transform plane, acting as a spatial frequency filter. Mathematically,

As a result of the spatial filter function F(2), the imaging resolution of Eq. 9 is much narrower than that of a first-order image; however, it is difficult to simplify this equation further in this form.

To get a better understanding of the physics behind Eq. 9, instead of modeling the finite radius of the lens as a disk, which results in the somb-function, we approximate the finite radius of the lens as the Gaussian function with a half-width D/2, and evaluate in one dimension. This leads to a Gaussian imaging-forming function instead of the somb-function. Working in one dimension, we change to x; 1 to x1; to x2, etc. Then ΔR(x1) simplifies to:

Corresponding to the experimental measurement, where D2 was placed at two off-center points in the Fourier transform plane, we model the filter function in one dimension by two rectangle functions: F(κ2) = Π(x − x2+, D) + Π(x − x2−, D), where , while A(x) = Π(x, w) + Π(x − 2w, w) + Π(x + 2w, w); assuming x2, x, and κ2(x) are integrated from ±∞, the resulting equation is an analytic expression. Define

which contains the imaging equation in Gaussian form for x1 and x2, in addition to some phase terms which are not observable in a first-order image, and

Then Eq. 10 is, after evaluating the κ2 and x2 integrals,

Then it is easy to see that restricting the allowed spatial frequencies of the Erfi functions constrains the values x is allowed to take, which, together with f(x1, x), improves the ability to resolve different points on the object plane. However, without evaluation Eq. 13 may still not be clear enough to show exactly how the resolution is affected, so we have included the following figures which plot some informative values to support our experimental observation. Figure 3(a) compares the theoretical first-order unresolved image of three slits with the second-order fluctuation correlation image calculated using Eq. 13. The filter function is calculated for a fiber diameter D of 0.105 mm and varying distances x from the center of the Fourier transform plane. Note that the plot for x = 0.18 mm demonstrates similar behavior to the observed experimental data, including the shift of the left and right peaks away from center. It is clear that, for a lens diameter of 1.36 mm, the gold transparent plot is completely unresolved. However, as the cutoff frequency κ2 = kx/f increases, the second-order resolution also increases as seen in the increasing separation of the peaks in Fig. 3. This is more clearly illustrated in Fig. 3(b) and (c); in (b) the imaging function of the first-order image (black) is plotted with the second-order imaging function (teal) for x = 0.18 mm. It is clear that the second-order imaging function is much narrower. In (c) the half-width of the imaging function in Eq. 13 is calculated at x1 = 0 for increasing values of x. Using the estimated experimental parameters, the second-order imaging resolution starts equal to the first order at x = 0 and increases to the experimental setup’s limit of about 0.002 mm.

Figure 3

All plots were calculated from Eq. 13.

(a) Theoretical comparison of the resolution of a 1.36 mm lens for three slits with a width and separation of 0.01241 mm. The transparent gold shading shows the first-order unresolved image, while the light gray indicates the ideal image; the lines plot Eq. 13 with spatial filter F(κ2) for several cutoff values , x2± = ±x, where k is the wavenumber and f is the focal length of lens L, versus the position of D1 in the image plane. As the cutoff frequency increases the ability to resolve the three slits also increases. (b) Comparison of the first-order Gaussian imaging function (black) to the second-order imaging function (teal) for x = 0.18 mm. The y-axis is in arbitrary units of normalized intensity. (c) Theoretical comparison of the resolution of Eq. 13 for a 1.36 mm lens as the cutoff frequency x increases. The dot marks the second-order imaging resolution of x = 0.18 mm from (a); the first-order imaging resolution corresponds to x = 0.

It is evident from the experiment and theoretical calculations that the increase in spatial resolution is strongly dependent on the chosen spatial filter. It is, in effect, applying a high-pass spatial filter to one arm, producing an edge-sharpening effect1112. The interesting part is that the correlation of the spatially filtered intensity fluctuations with arm one produces a resolved image, especially since neither arm “sees” a resolved image. This correlation filters out the lower spatial frequencies of the unresolved image of arm one, yielding a resolved image in the intensity fluctuation correlation rather than the intensity.

Conclusion

In summary, by using a high-pass spatial filter in the non-resolving side of a two-arm camera, the measurement of the intensity fluctuation correlation 〈ΔI1(1)ΔI2(2)〉 was able to resolve an object that could not be resolved by a traditional camera. This imaging method would be particularly useful for long-distance imaging in situations where it is impractical to have large lenses but high resolution is still desired, as it could take advantage of the large angular size of the sun, 0.5° relative to the earth, and the correspondingly small coherence length, on the order of 0.2 mm. In addition, since the thermal light image in 〈ΔI1(1)ΔI2(2)〉 is in general turbulence-free14, this method would be particularly attractive for satellite cameras taking high resolution images of objects on the ground. Technically more complicated optics or electronics for practical sunlight imaging will be discussed separately.

Additional Information

How to cite this article: Sprigg, J. et al. Super-resolution imaging using the spatial-frequency filtered intensity fluctuation correlation. Sci. Rep. 6, 38077; doi: 10.1038/srep38077 (2016).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.