Wigner function measurement using a lenslet array

Lei Tian; Zhengyun Zhang; Jonathan C. Petruccelli; George Barbastathis

doi:10.1364/OE.21.010511

1. Introduction

In order to extract the most information possible from an optical field, one can look for complete descriptions of the wave properties of a field that allow physically measurable quantities (e.g. intensity, energy density, Poynting vector) to be calculated over a region of free space. For most propagating optical fields, a description on some plane transverse to the optical axis is sufficient to determine the rest of the field; we will denote position on the plane by r = (x, y) and position along the optical axis by z. For monochromatic coherent light, the complex–valued field U(r; z) provides such a description. For a quasi–monochromatic partially coherent field, a common description is the mutual intensity, J(r₁, r₂; z), i.e. the statistical correlation of the field at all possible pairs of points over a plane [1]. The mutual intensity can be written as

J (r_{1}, r_{2}; z) = 〈 U (r_{1}; z) U^{*} (r_{2}; z) 〉,

where 〈·〉 denotes the ensemble average over statistical realizations of the field U, and the intensity is simply given by I(r) = J(r, r).

The mutual intensity may be measured directly by using interferometry: two sheared versions of a field are overlapped in a Young [2], Michelson stellar [3] or rotational shear [4,5] arrangement, and the fringe visibility gives the mutual intensity. However, interferometric methods are subject to mechanical scanning stability and signal–to–noise limitations. Other techniques use the fact that the intensity at any point in space may be computed from a pair of propagation integrals acting on J. By measuring the intensity over many different defocused planes, these integrals may be inverted to recover J through a process known as phase space tomography [6–11].

An alternative to the mutual intensity for describing an optical field is the radiance, a function that characterizes the distribution of power over position and direction. Let the radiance at some plane z be denoted by B(r, p; z) where p = (p_x, p_y) is the transverse component of the unit direction vector. Propagation of radiance is based on geometric optics and is simple to calculate:

B (r, p; z) = B (r - z p, p; 0)

The intensity at any point is given by the integral of B over all directions

I (r; z) = \int B (r, p; z) d^{2} p = \int B (r - z p, p; 0) d^{2} p .

Use of the radiance predates the wave theory of light, and it was initially described by assigning non-negative values to all trajectories coming from source points. Such descriptions are insufficient to model wave effects since these trajectories contain no information about constructive or destructive interference. However, in certain situations, these wave effects can be safely ignored, and a lenslet array can then be used to obtain an estimate for the radiance since they allow joint measurement of the spatial and directional distribution of light, as in Shack-Hartmann sensors [12], integral imaging systems [13–15] and light field cameras [16, 17].

More recently, generalized radiance functions have been proposed that describe interference while satisfying Eqs. (2) and (3); B is allowed to become negative to account for destructive interference, and it is also allowed to have non–zero values outside of the source [18]. The generalized radiance and mutual intensity both contain equivalent information: when defined over a plane, they are 4D functions that fully describe the intensity of an optical field (including wave effects) over a region of free space. In fact, for a quasi–monochromatic paraxial field of mean wavelength λ, the two descriptions are related by a simple Fourier transform relationship [19–21]:

\frac{1}{λ^{2}} B (r, λ u) = 𝒲 (r, u) = \int J (r + r^{'} / 2, r - r^{'} / 2) e^{- i 2 π u \cdot r^{'}} d^{2} r^{'},

where 𝒲 is known as the Wigner distribution function (WDF); it is related to the radiance by a change of variables from transverse direction vector p to spatial frequency u = (u_x, u_y) = p/λ. From Eq. (4), it is evident that if 𝒲 is known, an inverse Fourier transform can be used to obtain the mutual intensity.

In this manuscript we will show that not only can lenslet arrays measure radiance in the geometric optics sense, they can also provide information about the WDF when wave optics is considered. In order to understand how a lenslet array measurement relates to the WDF, we first consider a single lenslet centered at position r₀ and placed one focal length f away from a detector [22]. The effect of the lenslet’s aperture is to window a portion of the incident field centered at r₀. The intensity I(r) at the detector plane measures a spatial frequency distribution of the windowed field due to the Fourier transforming property of the lenslet. As illustrated in Appendix A, I (r) is related to the incident WDF 𝒲_i and the aperture WDF 𝒲_p by

I (r) = \frac{1}{λ f} S [𝒲_{i}, 𝒲_{p}] (r_{0}, u),

where spatial frequency u is mapped to position r at the detector plane according to u = (r − r₀)/λf, and S[𝒲₁, 𝒲₂](r, u) denotes a smoothed WDF (SWDF) defined over space and spatial frequency variables r and u, respectively, constructed from a convolution between the two WDFs 𝒲₁ and 𝒲₂ enclosed in the square brackets [22–26]:

S [𝒲_{1}, 𝒲_{2}] (r, u) = \int 𝒲_{1} (r^{'}, u^{'}) 𝒲_{2} (r^{'} - r, u - u^{'}) d^{2} r^{'} d^{2} u^{'} .

Equation 5 implies that for a lenslet centered at r₀ the detector pixel at r corresponds to the unique point at [r₀, (r − r₀)/λf] of the SWDF, as illustrated in Fig. 1(a). In order to fully measure the SWDF, one could scan the lenslet and take multiple images [24, 26].

Fig. 1 Illustration of the measurement under a single lenslet of the SWDF in one spatial dimension (r₀ = x₀ and u = u_x). (a) A lenslet centered at x₀ maps points (x₀, u_x) in the SWDF domain to positions x₀ + λfu_x at the lenslet’s Fourier plane. (b) A convolution of the aperture WDF with the incident WDF forms the SWDF.

Download Full Size | PDF

Since our goal is to accurately measure the incident WDF, we can interpret the construction of the SWDF in Eq. (5) as smoothing the incident WDF by convolution with the aperture WDF, as illustrated in Fig. 1(b). Ideally, then, one would seek to use an aperture for which 𝒲_p ∝ δ(r)δ(u), where δ is the Dirac delta function in order to measure the incident WDF with perfect resolution. However it may be easily verified from Eq. (4) that the existence of this aperture WDF is physically impossible; uncertainty relationships place limits on the product of the widths of 𝒲_p in space and spatial frequency. The finite extent of 𝒲_p limits the resolution in space and spatial frequency of the measurement of 𝒲_i. Despite these limitations, in many cases the SWDF itself is useful to provide a direct estimate of the WDF [24, 26]. However, if the exact WDF is required, Eq. (6) shows that the recovery of the WDF from the SWDF generally requires deconvolution [25].

A periodic array of lenslets enables measurements of different r₀ simultaneously and removes the need to scan. The advantage of using a lenslet array is that a detector pixel under a given lenslet at position r₀ measures a point (r₀, u) in the SWDF domain according to the mapping shown in Fig. 1(a). Therefore the geometry of an array of lenslets can be tailored to measure a desired region of the SWDF domain. In section 2 we rigorously derive the relationship between measured intensity and the SWDF of the incident field for an array of lenslets. For simplicity, we consider scalar fields in one spatial dimension. We demonstrate that, in general, the unique mapping implied by Eq. (5) no longer holds. We show that the intensity at a detector pixel in general contains light from multiple lenslets which we call cross–talk. Accurate measurement of the SWDF requires minimizing this cross–talk. In addition, both fully incoherent and fully coherent cases can have considerable amounts of cross–talk. In Section 3, we illustrate tradeoffs between coherence and fidelity using a numerical example, showing that there exists an optimal “Goldilocks” regime for array pitch, given the the coherence width of the input light, such that cross–talk is reduced to a minimum without the need for additional barriers to block light between lenslets. It is in this optimal regime that each detector pixel corresponds to a single point in the SWDF domain, allowing lenslet array systems to measure the SWDF with high accuracy. Since our goal is direct measurement of the field’s coherence properties, we also consider the discrepancy between the SWDF and WDF for this example.

2. Theory

We now present a rigorous analysis of a 1D field passing through a 1D lenslet array. This analysis can be easily extended to a 2D rectangular array; other configurations, such as a 2D hexagonal array, require straightforward modifications. Consider a quasi–monochromatic paraxial field with wavelength λ incident upon an ideal lenslet array with 100% fill factor, as illustrated in Fig. 2. A detector is placed at the back focal plane of the lenslet array, and we will refer to the region directly behind each lenslet on the detector as that lenslet’s detector cell. We assume an array of 2N + 1 identical unaberrated thin lenses, each of width w and focal length f. The transmittance function of such an array is given by

T (x) = \sum_{l = - N}^{N} rect (\frac{x - l w}{w}) exp [- \frac{i π}{λ f} {(x - l w)}^{2}],

where rect(·) denotes a rectangular function. In this configuration, we assign an integer index l to each lenslet, with the center lenslet having index l = 0; the N lenslets above and N lenslets below the center lenslet take on positive and negative values of l, respectively. Thus, the center of each lenslet is located at x = lw. We have assumed an odd number of lenslets to simplify notation, although the results we obtain can easily be extended to an even number of lenslets.

Fig. 2 Lenslet array geometry.

Download Full Size | PDF

The mutual intensity immediately to the right of the lenslet array is [1]

J_{1} (\bar{x} + \frac{x^{'}}{2}, \bar{x} - \frac{x^{'}}{2}) = J_{i} (\bar{x} + \frac{x^{'}}{2}, \bar{x} - \frac{x^{'}}{2}) T (\bar{x} + \frac{x^{'}}{2}) T^{*} (\bar{x} - \frac{x^{'}}{2}),

where x̄ and x′ are the center and difference coordinates, respectively, J_i is the mutual intensity of the illumination immediately before the lenslet array; the subscript i indicates that its associated function describes properties of the incident field at the input plane, and we will use this notation through the rest of the manuscript.

As a stepping stone to the full relationship between the incident field and the observed intensity behind the lenslet array, we will first consider a simpler system wherein we scan through the lenslets. That is, instead of letting light pass simultaneously through all the lenslets while recording the intensity image, we only let light pass through one lenslet at a time, cycling through all the lenslets while still recording a single image. This removes the effect of cross–lenslet interference, whose derivation we will consider later.

According to Eq. (5), each measurement samples the SWDF over spatial frequency with position fixed at the lenslet’s center, lw. The aperture of each lenslet is a rect function of width w, and thus the aperture WDF is given by

𝒲_{p} (x, u) = \frac{sin [2 π u (w - 2 | x |)]}{π u} rect (\frac{x}{w})

The total intensity at the detector plane is given by

I_{0} (x_{o}) = \frac{1}{λ f} \sum_{l = - N}^{N} S [𝒲_{i}, 𝒲_{p}] (l w, \frac{x_{o} - l w}{λ f}) .

where

\frac{1}{λ f} S [𝒲_{i}, 𝒲_{p}] (l w, \frac{x_{o} - l w}{λ f})

is the contribution of light through a single lenslet.

It is clear from this equation that the SWDF is sampled spatially at intervals of w, the spacing of the lenslet centers. The sampling rate along the spatial frequency axis in the SWDF is determined by both the detector pixel size and the linear mapping u = (x_o −lw)/λf between detector coordinate x_o and spatial frequency coordinate u. The mapping can be explained by the fact that (x_o − lw)/f equals to the angle between the ray reaching the detector pixel at x_o and the optical axis of the l^th lenslet under a small angle approximation. Note that if the angular spread of the SWDF is large enough, each detector cell will include contributions to intensity not only from the SWDF associated with its lenslet, but also from neighboring lenslets. This can be prevented by increasing the size of the lenslets or by decreasing the angular spread of the incident field by placing either a main lens with finite numerical aperture in front of the array [17] or physical barriers between lenslets [27]. If we assume that each detector cell measures only light from its associated lenslet, then we would have a detected intensity of the following form

I_{SWDF} (x_{o}) = \frac{1}{λ f} \sum_{l = - N}^{N} S [𝒲_{i}, 𝒲_{p}] (l w, \frac{x_{o} - l w}{λ f}) rect (\frac{x_{o} - l w}{w}) .

We refer to this expression as I_SWDF, because the intensity measured at x_o maps uniquely to the point [l̂w, (x_o − l̂w)/λf] in the SWDF, where l̂ is x_o/w rounded to the nearest integer.

The difference between Eq. (10) and Eq. (11) yields what we will call the 0^th order cross–talk term $I_{c}^{(0)}$ :

I_{c}^{(0)} (x_{o}) = \frac{1}{λ f} \sum_{l = - N}^{N} S [𝒲_{i}, 𝒲_{p}] (l w, \frac{x_{o} - l w}{λ f}) [1 - rect (\frac{x_{o} - l w}{w})] .

This term describes the additional light from other lenslets if the previous assumption were to not hold. In this situation multiple points in the SWDF can contribute to the same point on the detector.

To demonstrate the sampling described by Eqs. (10)–(12), an array containing three lenslets (centered at −w, 0, w) is shown in Fig. 3. According to Eq. (10), three lines sampled at spatial coordinates −w, 0, w parallel to the u–axis from the SWDF are mapped to the detector plane (marked by different colors in Fig. 3). To ensure one–to–one mapping, the maximum spatial frequency u_m of the l^th line sample cannot exceed w/2λf, as shown in case (a); otherwise, points at (lw, u_m) and [(l + 1)w, u_m − w/(λf)] from the SWDF domain will be measured by the same detector pixel at x_o = lw +λfu_m, as shown in case (b).

Fig. 3 Sampling of the SWDF using an array of three lenslets. (a) One–to–one mapping from the SWDF to the detector coordinate according to u = (x_o − lw)/(λf) as the angular spread of the SWDF is narrower than the numerical aperture of a lenslet. (b) Multiple points in the SWDF domain contribute to detector pixels in the cross–talk region as the angular spread of the incident field is wider than the numerical aperture of a lenslet, which produces the 0^th order cross–talk.

Download Full Size | PDF

So far, we have only considered the incoherent superposition of light from all of the lenslets, whereas light passing through all of the lenslets simultaneously should create additional interference terms. Since light from lenslets separated by a distance greater than the incident field’s coherence width will not create appreciable interference when mixed, it is useful to enumerate these cross–talk terms with an index n proportional to the lenslet separation. All possible pairs of lenslets with indices l′ and l″ such that |l′ − l″| = n > 0 contribute to the n^th order cross–talk term $I_{c}^{(n)}$ , given by

I_{c}^{(n)} (x_{o}) = \frac{2}{λ f} \sum_{l = - N + \frac{n}{2}}^{N - \frac{n}{2}} \iint 𝒲_{i} (\bar{x}, u) 𝒲_{p} (\bar{x} - l w, u - \frac{x_{o} - l w}{λ f}) cos [2 π (\frac{x - x_{o}}{λ f} + u) n w] d \bar{x} d u .

Note that when n is odd, l takes a value halfway between two integers, and thus 𝒲_p is centered at the edge between the (l −1/2)^th and (l +1/2)^th lenslets; when n is even, l takes every integer value, thus 𝒲_p is centered at the l^th lenslet. We expect the n = 1 term to be significant even in highly incoherent fields, since some points near the boundary between two neighboring lenslets are expected to be within the coherence width of the field.

The total output intensity, considering all of the discussed effects, can be written as the sum of three components

I (x_{o}) = I_{SWDF} (x_{o}) + I_{c}^{(0)} (x_{o}) + \sum_{n = 1}^{2 N} I_{c}^{(n)} (x_{o}) .

A detailed derivation of this result, obtained by performing Fresnel propagation integrals on Eq. (8), is given in Appendix B. Equation 14 demonstrates that if all orders of the cross–talk could be made small, then the measured intensity would be an accurate representation of the SWDF. In order for the cross–talk to be negligible, both the angular spread of the SWDF should be small [for

I_{c}^{(0)} (x_{o})

] and the coherence width should be less than the width of a single lenslet [for

I_{c}^{(n)} (x_{o})

]. In order to optimally measure the SWDF, the angular and coherence widths of the SWDF should be balanced so that as much of each lenslet’s detector cell is utilized as possible while minimizing cross–talk. It should also be noted that even with minimal cross–talk, the measurement yields only the SWDF; recovery of the mutual intensity (or WDF) of the field still requires deconvolution of the SWDF with the aperture WDF.

3. Numerical example

We study the effect of coherence width on the quality of the resulting measurement by the following example. Let us consider a spatially homogeneous incident field of wide enough spatial extent compared to the width of the lenslet array that it can be approximated as infinitely wide. The transverse coherence of the field is described using a Gaussian–correlated Schell model, such that the mutual intensity is

J_{i} (x_{1}, x_{2}) = exp [- \frac{{(x_{2} - x_{1})}^{2}}{2 σ_{c}^{2}}] .

The coherence width is proportional to the standard deviation σ_c of the coherence term. The WDF of the incident field is therefore independent of x̄ and is given by

𝒲_{i} (\bar{x}, u) = \frac{1}{\sqrt{2 π} σ_{u}} exp [- \frac{u^{2}}{2 σ_{u}^{2}}],

where σ_u = 1/(2πσ_c) quantifies the spatial frequency bandwidth of the WDF and is proportional to the angular spread of the field. The SWDF resulting from the convolution between the WDF of the input field and that of a rectangular aperture is

S [𝒲_{i}, 𝒲_{p}] (\bar{x}, u) = \frac{w^{2}}{\sqrt{2 π} σ_{u}} \int exp [- \frac{{(u - u^{'})}^{2}}{2 σ_{u}^{2}}] {[\frac{sin (π w u^{'})}{π w u^{'}}]}^{2} d u^{'} .

Therefore the SWDF term is

\begin{array}{l} I_{SWDF} (x_{o}) & = & \frac{w^{2}}{\sqrt{2 π} σ_{u} λ f} \sum_{l = - N}^{N} rect (\frac{x_{o} - l w}{w}) \int exp {- \frac{{[(x_{o} - l w) / λ f - u^{'}]}^{2}}{2 σ_{u}^{2}}} \\ [\frac{sin (π w u^{'})}{π w u^{'}}] d u^{'}, \end{array}

and the 0^th order cross–talk term is

\begin{array}{l} I_{c}^{(0)} (x_{o}) & = & \frac{w^{2}}{\sqrt{2 π} σ_{u} λ f} \sum_{l = - N}^{N} [1 - rect (\frac{x_{o} - l w}{w})] \\ \int exp {- \frac{{[(x_{o} - l w) / λ f - u^{'}]}^{2}}{2 σ_{u}^{2}}} {[\frac{sin (π w u^{'})}{π w u^{'}}]}^{2} d u^{'} . \end{array}

The n^th order cross–talk term by carrying out the integration in Eq. (13) is

\begin{array}{l} I_{c}^{(n)} (x_{o}) & = & \frac{2 w^{2}}{\sqrt{2 π} σ_{u} λ f} \sum_{l = - N + \frac{n}{2}}^{N - \frac{n}{2}} \int exp {- \frac{{[(x_{o} - l w) / λ f - u^{'}]}^{2}}{2 σ_{u}^{2}}} cos (2 π n w u^{'}) \\ \frac{sin [π w (n u_{0} + 2 u^{'}) / 2]}{π w (n u_{0} + 2 u^{'}) / 2} \frac{sin [π w (n u_{0} - 2 u^{'}) / 2]}{π w (n u_{0} - 2 u^{'}) / 2} d u^{'}, \end{array}

where u₀ = w/λf is the spatial frequency support of a single lenslet. The total output intensity measured at a detector pixel will be the sum of the SWDF and all cross–talk terms, as given in Eq. (14).

Three different cases with varying coherence widths are simulated based on the results in Eqs. (17)–(20). In the simulation, the wavelength of the incident field is 500nm. Five lenslets are used in the array, each having width w = 330μm and focal length f = 5mm, yielding spatial frequency support of u₀ = 0.132μm⁻¹. The simulation results are shown in Fig. 4. For all three cases, the total output intensity in row (a) is composed of the SWDF term in row (b) and the total contribution of cross–talk in row (c). The total cross–talk is further analyzed by decomposing it as the 0^th order term in row (d) and the total of higher order terms in row (e). Simulations on arrays with larger numbers of lenslets were also conducted; results are not shown here because they are very similar to the ones in Fig. 4.

Fig. 4 Left: highly incoherent; middle: highly coherent; and right: partially coherent case. (a) Total output intensity is composed of (b) SWDF term and (c) total contribution from cross–talk terms. The total cross–talk is composed of (d) 0^th order cross–talk and (e) total of higher order cross–talk. All the intensities are normalized to the maximum value in the total output. The horizontal axis is the spatial coordinate normalized by the width of a lenslet.

Download Full Size | PDF

In the highly incoherent case shown in the left column (σ_c = 0.01w), higher order cross–talk is minimal. However, due to the large angular spread in the incident field, the measurement is corrupted by 0^th order cross–talk. The opposite is the highly coherent case, shown in the middle column (σ_c = 20w). Here, most of the cross–talk comes from higher order terms. The results for a partially coherent field (σ_c = 0.1w) is shown in the right column; cross–talk contributes minimally to the final intensity, although both 0^th order and higher order terms are present.

Notice that although the SWDF itself is homogeneous in x, the intensities measured in Fig. 4 are not identical under each lenslet. The reason for this is cross–talk. As expressed in Eq. (14) cross–talk under a given lenslet results from the light from neighboring lenslets. The number of neighboring lenslets contributing to the cross–talk depends on both the coherence width of the field and on the angular spread of the field. For a lenslet near the edge of the array, there are few neighbors in the direction of the edge, meaning less cross–talk from that direction. Since our simulated lenslet array consists of only 5 lenslets for simplicity, at least 2 of the lenslets are “edge lenslets,” and depending on the coherence width and angular spread, edge effects may influence all lenslets. Practical lenslet arrays are likely to contain considerably more lenslets, and generally most lenslets will have minimal contributions from edge effects. Because of this, in comparing measured intensity to the WDF and SWDF of the incident field, we consider only the intensity under the central lenslet in our simulated array.

The spatial position under the central lenslet maps into the SWDF domain according to (x = 0, u = x_o/λf). In Fig. 5 we compare the total output intensity under the central lenslet (dotted green lines) to the corresponding slice (u = x_o/λf) of the SWDF (dashed blue lines) for each of the three simulated fields. The matching slice of the WDF (solid red lines) illustrates the effect of aperture convolution that generates the SWDF. In both the highly incoherent and partially coherent cases, the SWDF and WDF are very similar, since the WDF of the aperture is much smaller than any variations in the incident WDFs. In the highly coherent case, the incident WDF is narrower in u than the aperture WDF, and therefore the SWDF is significantly broadened by the convolution. In order to recover the WDF from the measured intensity, deconvolution is necessary [25].

Fig. 5 Comparison of WDF (solid red line), SWDF (dashed blue lines) and total output intensity (dotted green lines) for (a) highly incoherent (σ_c = 0.01w), (b) highly coherent (σ_c = 20w), and (c) partially coherent (σ_c = 0.1w) incident light.

Download Full Size | PDF

We quantify how close the lenslet measurement is to the WDF of the incident field by comparing the total output intensity under the central lenslet to the incident WDF at those (x, u) coordinates using an error metric R_error:

R_{error} = \frac{{| total output intensity - incident WDF |}^{2}}{{| total output intensity |}^{2}} .

The total output intensity contains contributions both from the SWDF and from cross–talk. To compare the relative importance of these contributions, we quantify the cross–talk corruption in the output through the cross–talk intensity fraction R_cross–talk as

R_{cross - talk} = \frac{{| total intensity in all terms of cross - talk |}^{2}}{{| total output intensity |}^{2}} .

The SWDF itself is a smoothed version of the incident WDF. To quantify the discrepancy due to this smoothing, we define the convolution error R_conv as

R_{conv} = \frac{{| SWDF - incident WDF |}^{2}}{{| total output intensity |}^{2}} .

All these error metrics are plotted as functions of the coherence width of the incident field normalized to the lenslet width (σ_c/w) in Fig. 6. As seen in the dashed green curve, the contribution from cross–talk increases quickly as the field becomes less coherent. When the field becomes more coherent, the contribution from cross–talk also increases until it saturates to the point in which the field is coherent within the whole array. There exists a partially coherent regime where the SWDF can be measured with minimal cross–talk corruption. Depending on accuracy requirements, this regime may provide acceptable measurements. For example, if less than 1% of cross–talk can be tolerated, then the coherence width should be such that 0.02w < σ_c < w. On the other hand, the convolution error increases as the field becomes more coherent, making the SWDF a less accurate estimate of the WDF in these situations. R_error considers artifacts from both cross–talk and convolution, has a similar shape to the cross–talk curve. The measurement deviates from the original WDF except in a partially coherent region. If error needs to be at most 1%, then we would need 0.02w < σ_c < 0.4w.

Fig. 6 R_error in solid blue curve, cross–talk intensity fraction R_cross–talk in dashed green curve, and convolution error R_conv in red dotted curve as functions of the normalized coherence length of incident light σ_c/w.

Download Full Size | PDF

4. Concluding Remarks

Although the numerical example was chosen explicitly to consider the effect of coherence width on the measurement of the SWDF using a lenslet array, this simple model can also provide useful insights for a much broader class of fields whose intensity varies slowly across the field, with features much wider than the coherence width. As a rule of thumb, higher order (coherent) cross–talk can be reduced by ensuring that the lenslet apertures are at least one coherence width in size. This makes intuitive sense, since an aperture larger than the coherence width will not cause the incident beam to diffract significantly, and any light that is diffracted from the aperture will not interfere with that from neighboring lenslets. Both 0^th and higher order cross–talk can be reduced by ensuring the incident illumination’s angular spread is such that each lenslet primarily illuminates only the pixels lying within its detector cell, such that there is a nearly one–to–one mapping from SWDF space to each detector pixel.

It should also be noted that we have derived these results under the paraxial approximation and that both the 0^th and higher order cross–talk can include contributions for which the light propagates highly non–paraxially from one lenslet to its neighbors. In these cases, we expect that a similar analysis can be performed using non–paraxial versions of the Wigner function [28, 29], although this is outside the scope of our current work.

As was discussed while analyzing the example, there are cases where the SWDF is not an accurate estimate of the WDF. As illustrated in Fig. 1(b), the SWDF can be interpreted as blurring the incident field’s WDF by convolution with the aperture WDF. The problem of recovering the incident WDF is therefore similar to improving the resolution of an optical imaging system through deconvolution. We therefore expect there to be significant difficulties especially when the incident WDF contains features smaller than the aperture WDF. Performing deconvolution to recover the WDF may benefit from techniques such as coded apertures [30, 31] and compressed sensing [10, 32].

Appendix A: Derivation of Eqs. (5) and (6)

Consider an incident field with mutual intensity J_i passing through a lenslet centered at r₀ whose aperture function is given by p. The mutual intensity immediately after the lenslet is

J (r + \frac{r^{'}}{2}, r - \frac{r^{'}}{2}) = J_{i} (r + \frac{r^{'}}{2}, r - \frac{r^{'}}{2}) T (r + \frac{r^{'}}{2}) T^{*} (r - \frac{r^{'}}{2}),

where T(r) = p(r − r₀) exp[−iπ|r − r₀|²/(λf)] is the transmittance function of the lenslet. It is straightforward to verify that the WDF of the product of two mutual intensities is equivalent to the convolution along u of the WDFs of the two mutual intensities. Using this property of WDFs, we can express the WDF immediately after the lenslet 𝒲_o as:

𝒲_{o} (r, u) = \int 𝒲_{i} (r, u^{'}) 𝒲_{T} (r, u - u^{'}) d^{2} u^{'},

where 𝒲_T is the WDF of the transmittance function and all WDFs described here are in the lenslet plane. Propagation to the detector plane at z = f is performed using Eq. (2), such that the detector plane WDF 𝒲_d is given by

𝒲_{d} (r, u) = 𝒲_{o} (r - λ f u, u) = \int 𝒲_{i} (r - λ f u, u^{'}) 𝒲_{T} (r - λ f u, u - u^{'}) d^{2} u^{'} .

By using the definition of T and Eq. (4), 𝒲_T is related to the aperture WDF 𝒲_p by

𝒲_{T} (r, u) = 𝒲_{p} (r - r_{0}, \frac{r - r_{0}}{λ f} + u),

which can be substituted into Eq. (26) to yield

𝒲_{d} (r, u) = \int 𝒲_{i} (r - λ f u, u^{'}) 𝒲_{p} (r - λ f u - r_{0}, \frac{r - r_{0}}{λ f} - u^{'}) d^{2} u^{'} .

The intensity at the detector plane is given by an integral of 𝒲_d over u [Eq. (3)], yielding

I (r) = \int 𝒲_{i} (r - λ f u, u^{'}) 𝒲_{p} (r - λ f u - r_{0}, \frac{r - r_{0}}{λ f} - u^{'}) d^{2} u^{'} d^{2} u,

which proves Eqs. (5) and (6).

Appendix B: Proof of Eq. (14)

The intensity at the output plane x_o at one focal length to the right of the lenslet array is related to J₁ by double Fresnel integrals under the paraxial approximation,

I (x_{o}) = \frac{1}{λ f} \iint J_{1} (x_{1}, x_{2}) exp {\frac{i π}{λ f} [{(x_{o} - x_{1})}^{2} - {(x_{o} - x_{2})}^{2}]} d x_{1} d x_{2} .

Substitutions of Eqs. (7) and (8) into Eq. (30) shows that the calculation of I(x_o) requires in general carrying out double summations with respect to different lenslet indices l₁ and l₂. After some simplification, the output intensity is rewritten as

\begin{array}{l} I (x_{o}) = \frac{1}{λ f} \sum_{l_{1} = - N}^{N} \sum_{l_{2} = - N}^{N} \iint J_{i} (x_{1}, x_{2}) rect (\frac{x_{1} - l_{1} w}{w}) rect (\frac{x_{2} - l_{2} w}{w}) \\ exp {- i \frac{2 π}{λ f} [(l_{1}^{2} - l_{2}^{2}) \frac{w^{2}}{2} + (x_{1} - x_{2}) x_{o} - l_{1} w x_{1} + l_{2} w x_{2}]} d x_{1} d x_{2} . \end{array}

The double summations are simply due to the fact that the intensity at the region directly behind a particular lenslet results from both the field passing through the local lenslet and cross–talk from the field passing through neighboring lenslets. We proceed with the derivation by the change of variables

\bar{x} = \frac{x_{1} + x_{2}}{2}, x^{'} = x_{1} - x_{2},

and

m = l_{1} + l_{2}, n = l_{1} - l_{2} .

The consequence of the change of variables in Eq. (33) can be studied as follows. a) When n is even, we can write

n = 2 q, where q = - N, - N + 1, \dots, N - 1, N;

and since n = l₁ − l₂, we can assume

l_{1} = l + q, l_{2} = l - q,

and

m = 2 l, where l = - N + | q |, - N + | q | + 1, \dots, N - | q | - 1, N - | q | .

b) When n is odd, we can again write

n = 2 q, where q = - N + \frac{1}{2}, - N + \frac{3}{2}, \dots, N - \frac{3}{2}, N - \frac{1}{2},

and

l_{1} = l + q, l_{2} = l - q,

with

m = 2 l, where l = - N + | q |, - N + | q | + 1, \dots, N - | q | - 1, N - | q | .

The substitution of the change of variables leads

\begin{array}{l} I (x_{o}) & = & \frac{1}{λ f} \sum_{\begin{array}{l} n = - 2 N \\ even \end{array}}^{2 N} \sum_{l = - N + \frac{| n |}{2}}^{N - \frac{| n |}{2}} \iint rect (\frac{\bar{x} - l w + (x^{'} - n w) / 2}{w}) rect (\frac{\bar{x} - l w + (x^{'} - n w) / 2}{w}) \\ J_{i} (\bar{x} + \frac{x^{'}}{2}, \bar{x} - \frac{x^{'}}{2}) exp [- i \frac{2 π}{λ f} (x_{o} - l w) x^{'} + i \frac{2 π}{λ f} (\bar{x} - l w) n w] d \bar{x} d x^{'} \\ + & \frac{1}{λ f} \sum_{\begin{array}{l} n = - 2 N + 1 \\ odd \end{array}}^{2 N - 1} \sum_{l = - N + \frac{| n |}{2}}^{N - \frac{| n |}{2}} \iint rect (\frac{\bar{x} - l w + (x^{'} - n w) / 2}{w}) rect (\frac{\bar{x} - l w - (x^{'} - n w) / 2}{w}) \\ J_{i} (\bar{x} + \frac{x^{'}}{2}, \bar{x} - \frac{x^{'}}{2}) exp [- i \frac{2 π}{λ f} (x_{o} - l w) x^{'} + i \frac{2 π}{λ f} (\bar{x} - l w) n w] d \bar{x} d x^{'} . \end{array}

Notice that a term of fixed n contributes a non–zero value to I(x_o) only if the two rect–functions overlap. This implies that the separation x′ between the pair of correlating points on the incident field can only take certain values, as determined by the following inequalities

| \bar{x} - l w | < w / 4,

(n - 1) w + 2 | \bar{x} - l w | < x^{'} < (n + 1) w - 2 | \bar{x} - l w | .

Eq. (42) implies that x′ is bounded to a region of width 2w − 4|x̄ − lw| centered at nw. Also recall that the magnitude of mutual intensity is significantly larger than zero at large separation distance x′ only if the field is highly coherent. This implies that more terms in the summation over n need to be considered if the field is more coherent. To simplify Eq. (40), we relate I(x_o) to the WDF of the incident field and the WDF 𝒲_p(x̄, u) of a rectangular aperture of width w, by completing the integration with respect to x′ to yield

\begin{array}{l} I (x_{o}) = \frac{1}{λ f} \\ {\sum_{\begin{matrix} n = - 2 N + 1 \\ odd \end{matrix}}^{2 N - 1} \sum_{l = - N + \frac{| n |}{2}}^{N - \frac{| n |}{2}} \iint 𝒲_{i} (\bar{x}, u) 𝒲_{p} (\bar{x} - l w, u - \frac{x_{o} - l w}{λ f}) exp [i 2 π (\frac{\bar{x} - x_{o}}{λ f} + u) n w] d \bar{x} d u \\ + \sum_{\begin{matrix} n = - 2 N \\ even \end{matrix}}^{2 N} \sum_{l = - N + \frac{| n |}{2}}^{N - \frac{| n |}{2}} \iint 𝒲_{i} (\bar{x}, u) 𝒲_{p} (\bar{x}, l w, u - \frac{x_{o} - l w}{λ f}) exp [i 2 π (\frac{\bar{x} - x_{o}}{λ f} + u) n w] d \bar{x} d u .} \end{array}

Finally, by combining the complex conjugate terms in n, we arrive at Eq. (14).

Acknowledgments

The authors thank Laura A. Waller and Hanhong Gao for their helpful discussions. This research was funded by the Chevron Energy Technology Company and Foxconn Technology Group.

References and links

1. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995).

2. B. J. Thompson and E. Wolf, “Two-beam interference with partially coherent light,” J. Opt. Soc. Am. 47, 895 (1957) [CrossRef] .

3. W. Tango and R. Twiss, “Michelson stellar interferometry,” Prog. Optics 17, 239–277 (1980) [CrossRef] .

4. K. Itoh and Y. Ohtsuka, “Fourier-transform spectral imaging: retrieval of source information from three-dimensional spatial coherence,” J. Opt. Soc. Am. A 3, 94–100 (1986) [CrossRef] .

5. D. L. Marks, R. A. Stack, and D. J. Brady, “Three-dimensional coherence imaging in the Fresnel domain,” Appl. Opt. 38, 1332–1342 (1999) [CrossRef] .

6. M. G. Raymer, M. Beck, and D. McAlister, “Complex wave-field reconstruction using phase-space tomography,” Phys. Rev. Lett. 72, 1137–1140 (1994) [CrossRef] .

7. K. G. Larkin and C. J. R. Sheppard, “Direct method for phase retrieval from the intensity of cylindrical wave fronts,” J. Opt. Soc. Am. A 16, 1838–1844 (1999) [CrossRef] .

8. D. M. Marks, R. A. Stack, and D. J. Brady, “Astigmatic coherence sensor for digital imaging,” Opt. Lett. 25, 1726–1728 (2000) [CrossRef] .

9. S. Cho and M. A. Alonso, “Ambiguity function and phase-space tomography for nonparaxial fields,” J. Opt. Soc. Am. A 28, 897–902 (2011) [CrossRef] .

10. L. Tian, J. Lee, S. B. Oh, and G. Barbastathis, “Experimental compressive phase space tomography,” Opt. Express 20, 8296–8308 (2012) [CrossRef] [PubMed] .

11. L. Tian, S. Rehman, and G. Barbastathis, “Experimental 4D compressive phase space tomography,” in “Frontiers in Optics,” (Optical Society of America, 2012), p. FM4C.4.

12. B. C. Platt and R. Shack, “History and principles of Shack-Hartmann wavefront sensing,” J. Refract. Surg. 17(2001) [PubMed] .

13. G. Lippmann, “La photographie integrale,” Comptes-Rendus, Academie des Sciences 146 (1908).

14. A. Stern and B. Javidi, “Three-dimensional image sensing and reconstruction with time-division multiplexed computational integral imaging,” Appl. Opt. 42, 7036–7042 (2003) [CrossRef] [PubMed] .

15. J.-H. Park, K. Hong, and B. Lee, “Recent progress in three-dimensional information processing based on integral imaging,” Appl. Opt. 48, H77–H94 (2009) [CrossRef] [PubMed] .

16. E. H. Adelson and J. Y. A. Wang, “Single lens stereo with a plenoptic camera,” IEEE Trans. Pattern Anal. Mach. Intell. 14, 99–106 (1992) [CrossRef] .

17. R. Ng, M. Levoy, M. Bredif, G. Duval, M. Horowitz, and P. Hanrahan, “Light field photography with a hand-held plenoptic camera,” Tech. Rep. CTSR 2005-02, Stanford (2005).

18. A. T. Friberg, “On the existence of a radiance function for finite planar sources of arbitrary states of coherence,” J. Opt. Soc. Am. 69, 192–198 (1979) [CrossRef] .

19. E. Wigner, “On the quantum correction for thermodynamic equilibrium,” Phys. Rev. 40, 0749–0759 (1932) [CrossRef] .

20. L. Dolin, “Beam description of weakly-inhomogeneous wave fields,” Izv. Vyssh. Uchebn. Zaved. Radiofiz 7, 559–563 (1964).

21. A. Walther, “Radiometry and coherence,” J. Opt. Soc. Am. 58, 1256–1259 (1968) [CrossRef] .

22. Z. Zhang and M. Levoy, “Wigner distributions and how they relate to the light field,” in “IEEE International Conference on Computational Photography (ICCP),” (IEEE, 2009), pp. 1–10 [CrossRef] .

23. H. Bartelt, K.-H. Brenner, and A. Lohmann, “The Wigner distribution function and its optical production,” Opt. Commun. 32, 32–38 (1980) [CrossRef] .

24. A. Wax and J. E. Thomas, “Optical heterodyne imaging and Wigner phase space distributions,” Opt. Lett. 21, 1427–1429 (1996) [CrossRef] [PubMed] .

25. H. N. Chapman, “Phase-retrieval X-ray microscopy by Wigner–distribution deconvolution,” Ultramicroscopy 66, 153–172 (1996) [CrossRef] .

26. L. Waller, G. Situ, and J. Fleischer, “Phase–space measurement and coherence synthesis of optical beams,” Nat. Photonics 6, 474–479 (2012) [CrossRef] .

27. H. Choi, S.-W. Min, S. Jung, J.-H. Park, and B. Lee, “Multiple-viewing-zone integral imaging using a dynamic barrier array for three-dimensional displays,” Opt. Express 11, 927–932 (2003) [CrossRef] [PubMed] .

28. K. B. Wolf, M. A. Alonso, and G. W. Forbes, “Wigner functions for Helmholtz wave fields,” J. Opt. Soc. Am. A 16, 2476–2487 (1999) [CrossRef] .

29. S. Cho, J. Petruccelli, and M. Alonso, “Wigner functions for paraxial and nonparaxial fields,” J. Mod. Optic. 56, 1843–1852 (2009) [CrossRef] .

30. N. Lindlein, J. Pfund, and J. Schwider, “Algorithm for expanding the dynamic range of a shack-hartmann sensor by using a spatial light modulator array,” Opt. Eng. 40, 837–840 (2001) [CrossRef] .

31. M. E. Gehm, S. T. McCain, N. P. Pitsianis, D. J. Brady, P. Potuluri, and M. E. Sullivan, “Static two-dimensional aperture coding for multimodal, multiplex spectroscopy,” Appl. Opt. 45, 2965–2974 (2006) [CrossRef] [PubMed] .

32. Z. Zhang, Z. Chen, S. Rehman, and G. Barbastathis, “Factored form descent: a practical algorithm for coherence retrieval,” Opt. Express 21, 5759–5780 (2013) [CrossRef] [PubMed] .

Wigner function measurement using a lenslet array

Abstract

1. Introduction

2. Theory

3. Numerical example

4. Concluding Remarks

Appendix A: Derivation of Eqs. (5) and (6)

Appendix B: Proof of Eq. (14)

Acknowledgments

References and links

Cited By

Figures (6)

Equations (43)

Optics Express