Abstract
Under the scalar paraxial approximation, an optical wavefield is considered to be complex function dependent on position; i.e., at a given location in space the optical field is a complex value with an intensity and phase. The optical wavefield propagates through space and can be modeled using the Fresnel transform. Lenses, apertures, and other optical elements can be used to control and manipulate the wavefield and to perform different types of signal processing operations. Often these optical systems are described theoretically in terms of linear systems theory leading to a commonly used Fourier optics framework. This is the theoretical framework that we will assume in this manuscript. The problem which we consider is how to recover the phase of an optical wavefield over a plane in space. While today it is relatively straightforward to measure the intensity of the optical wavefield over a plane using CMOS or CCD sensors, recovering the phase information is more complicated. Here we specifically examine a variant of the problem of phase retrieval using two intensity measurements. The intensity of the optical wavefield is recorded in both the image plane and the Fourier plane. To make the analysis simpler, we make a series of important theoretical assumptions and describe how in principle the phase information can be recovered. Then, a deterministic but iterative algorithm is derived and we examine the characteristics and properties of this algorithm. Finally, we examine some of the theoretical assumptions we have made and how valid these assumptions are in practice. We then conclude with a brief discussion of the results.
1. Introduction
Visible light is an electromagnetic field with a wavelength that ranges roughly from ultraviolet (300 nm) to infrared (800 nm). An electromagnetic field is vectorial in nature with three electric and magnetic orthogonal field components. When designing a microwave radio receiver (long wave wavelengths relative to visible light) one must consider how these different magnetic and electric field components interact with each other and with conducting strips of metal and dielectric substrates that are used for impedance matching, mixing, and waveguiding. The resulting equations (that follow from Maxwell’s equations) are quite complex and require intensive numerical calculation procedures since the different vectorial components of the EM field interact with each other forming coupled vectorial equations [1].
Fortunately when examining optical problems in the visible regime it is often possible to make the scalar approximation. This means that we can assume that the different vectorial components of the optical wavefield do not “see” each other and hence each vectorial component can be treated independently as a separate scalar problem. This approximation is very useful and quite accurate and can be used to analyze a wide range of important and practical optical problems.
Therefore in this manuscript we assume that a scalar description of light propagation is valid and use the Fresnel transform to relate an optical field in one plane to that in another plane, where the optical planes are separated from each other by an axial distance . The light sources used are assumed to be of a definite monochromatic temporal frequency and are both temporally and spatially coherent [2, 3]. Optical elements used to shape the form of the illuminating wavefield are assumed to be “thin” and operate on an incident optical field, multiplying it by a function describing the optical element over a plane [4]. The assumptions we have described place us in a Fourier optics framework [4]. Fourier optics is an important branch of optical theory since it allows the development of relatively simple and intuitive models of optical systems. These models provide significant insight into the characteristics of the underlying optical systems and are also reasonably accurate. Furthermore Fourier optics builds important bridges between disciplines in particular allowing optical systems to be interpreted in terms of signal processing operations that are commonly used in electrical engineering, communication theory, and control engineering [4–9]. This type of viewpoint can be extended still further with the development of mixed space/time transforms or mixed space/spatial frequency transforms such as the fractional Fourier transform [10], the Linear Canonical Transform (LCT), and perhaps even more generally the form of Wigner distribution functions [11–17].
We note that in this theoretical framework if we know the intensity and phase of the optical field in a give optical plane, it is then possible to calculate (numerically or under certain circumstances analytically) the field distribution in any other plane using the LCT [18, 19] and hence we can image over a 3D volume, with varying levels of spatial resolution [19].
In a practical sense then it is important to be able to measure experimentally both the intensity and the phase of an optical wavefield in a given plane. Nowadays high quality low noise CCD and CMOS sensors are commercially available and widely used. With these devices it is relatively straightforward to measure the intensity of an optical field in a particular plane. We will now briefly discuss several techniques for indirectly measuring the phase distribution of the field.
1.1. Digital Holography
Digital holography (DH) was first proposed by Gabor himself in [20–22] and the first experimental implementation was done by Goodman et al. in [23]; later pioneering work was done by Yaroslawski, [24], Onural and Scott [25], and Schnars and Jüptner [26]. More recent theoretical and experimental analysis of these systems is given [19, 27–36]. In holography we need both an object wave and a so-called reference wave. Light from the laser source in the optical setup is divided by a beam splitter into a reference and illuminating plane waves. The illuminating plane wave is directed by a mirror so that it illuminates the object we wish to examine. The wavefront interacts with a “thin” object and the resulting object wavefield scatters and propagates a distance to the CCD/CMOS plane. Meanwhile the reference wave is directed through the optical system where it is recombined with the object wave at a second beam splitter. Thus at the sensor plane we record the interference between the reference and object wavefields thereby recording a hologram. A drawback of this approach is that the additional optics parts are required; a clean stable reference beam is essential; and several undesired “noise” terms, two DC terms and a twin image term, are produced by the recording process. These terms can be removed by capturing several phase shifted intensity images and performing some postprocessing operations. Single shot off-axis holography with a tilted plane reference wave is possible; however the spatial resolution of the resulting real image term (object wavefield) is significantly reduced. The accuracy and limitations of this approach have already been examined in several publications; see, for example, [33, 37].
1.2. Transport of Intensity
These approaches are quite different in nature and require the capture of two separate intensity distributions. Usually these distributions are captured in an image plane and in a slightly defocused plane. Using the so-called transport of intensity equation, see, for example, [38, 39], it is possible to recover an estimate for the phase distribution. There are several interesting aspects of the reconstruction of the phase; firstly the reconstructed phase is not wrapped over 2 and secondly with some reconstruction approaches, namely, an FFT based approach to indefinite integration, see [40], low frequencies can be more difficult to recover.
1.3. Iterative Phase Retrieval
Finally we turn our attention to iterative phase retrieval approaches. Here usually two or more intensity distributions are captured in different LCT domains, or in physical optics terms the intensity is captured in several different optical planes. Often these different optical planes are separated from each other by sections of free-space and combinations of lenses, but also in the image and corresponding Fourier transform planes. The origin of these types of algorithms dates back to the iterative Gerschberg-Saxton (GS) approach first published in the 70s and was further developed Fienup et al. [41, 42]. The initial algorithms assumed that two measurements were made and the intensity distribution was recorded in the image plane and in the Fourier plane. Later algorithms extended these approaches to work with multiple different Fresnel domains [43–47].
The image-Fourier plane (GS) algorithm works as follows:(1)One starts with an initial guess of the phase distribution in the image plane(2)We then numerically calculate the Fourier transform of this distribution and compare the calculated intensity in the Fourier plane to the experimentally measured result. The calculated intensity distribution will differ from the experimentally measured distribution because our initial guess at the phase distribution will not be correct(3)At this stage we replace the incorrect numerically calculated intensity distribution with the correct experimentally measured intensity, while retaining the phase values and numerically propagate the distribution back to the image plane(4)We do not expect the calculated image intensity to be the same as the measured intensity distribution. Again we replace the numerically calculated distribution with the measured distribution(5)Steps (2)-(4) are repeated until the algorithm converges, i.e., when there is little difference between the measured and calculated intensity distributions. The similarity between the measured and numerically calculated distributions can be estimated with a root mean square error measure [46].
In [41], Fienup directly compares the behavior of this algorithm to gradient descent optimization approaches. The behavior and properties of these algorithms have been of significant interest in the field. It is generally said that for phase retrieval to work effectively one needs to have a good estimate of the initial phase. When that is the case the algorithms tend to converge to the nearest local minimum which in this case would be one of the global minimums. And hence the algorithm will return the real-physical phase distribution and not a solution that merely satisfies the numerical constraints of the problem but which has no physical meaning. In this manuscript we will explore this question of multiple phase solutions more closely, by examining a simplified version of the phase retrieval problem. A major outstanding issue with these iterative approaches is that algorithm often stagnates in a local minimum yielding a poor estimate for the phase. There is also the possibility that there are several different solutions to this ill-posed problem, at least in the 1D case [48].
In the iterative phase retrieval approach that we have just been discussing the intensity distributions are captured using CCD/CMOS arrays. And so in practice the intensity distribution is only recorded over a finite spatial extent and with a finite number of pixels which limits the spatial resolution and also introduces numerical sampling effects [37, 49]. The numerical procedures that are then used with these 2D arrays of intensity values tend to rely on fast Fourier transform (FFT) based algorithms for their implementation. This requires developing the software and numerical techniques for implementing these algorithms so that they are sufficiently fast and so that other sampling rules and numerical issues can be addressed [49]. This development, while necessary, tends to concentrate attention on the software implementation rather than on other physical and mathematical modeling aspects of the problem.
Therefore in this manuscript we deliberately take a different approach that emphasizes a theoretical analysis using a relatively simple model. Similar to the iterative phase retrieval approach we analyze how the phase may be estimated from two intensity measurements: one made in the image plane and one made in the Fourier plane. We do however make an important simplification. We assume that we can measure the intensity distribution in the Fourier plane over an infinite extent and ignore any sampling or discretization process that occurs due to the finite extent of the pixels in the CCD/CMOS array. Furthermore we assume that the intensity distribution in the image plane consists of a finite set of discrete point sources. With these two simplifications we can derive an analytical solution to the phase retrieval problem.
The manuscript is organized in the following manner: In Section 2, we show, under assumptions, which we have outlined in the text above and which are detailed in Section 2, how the phase values can in principle be determined using an analytical solution. To recover the phase values we must use an iterative procedure which has a finite but large number of steps. In Section 3, we examine an algorithm for implementing the iterative strategy derived in Section 2 and examine how the number of operations depends on the number of contributing point sources in the image plane. We show that there are multiple solutions to the phase retrieval problem. We then examine approaches for reducing the number of possible solutions. In Section 4, we finish with a brief discussion of the presented analysis and highlight some practical difficulties when extending the theoretical model to realistic real-world phase measurements.
2. Analysis of an Idealized Experiment
We begin our analysis by examining the idealized optical system depicted in Figure 1. There is a set of distinct point sources, each with a known intensity value and hence a known magnitude, , where is the magnitude of point source . The phase value, , associated with each point source however remains the unknown physical parameter that we wish to measure. So to reiterate, the problem as depicted in Figure 1 is that we are given a set of point sources with known amplitudes and unknown phase values. We further assume that these points sources are uniformly spaced a distance from each other. We are also given a second measurement, the intensity distribution in the Fourier plane. What we have just described is very similar to the practical iterative phase retrieval experiment described in the Introduction. In a real experiment however these intensity distributions are recorded by CCD/CMOS arrays which have a finite extent, a finite number of pixels, and which sample of the intensity distribution they measure. These are critical factors when considering the performance limitations of these systems; see [33, 36]. They also complicate the theoretical analysis and so now we neglect these factors so that we can develop a simpler theoretical description of the phase retrieval problem.
We imagine that we can measure the intensity in the Fourier domain over an infinite extent and with a sampling step size that is extremely small, so that any sampling effects can be ignored. We also imagine that only a finite number of perfect point sources in the input plane (with known intensity but unknown phase) contribute to the Fourier plane intensity distribution. Hence in this analysis we exclude any physical imposed limiting factors on the ability of our system to make a measurement. This will allow us to significantly simplify the analysis and to concentrate on another and arguable more important effect: the existence of a very high number of alternative solutions to the phase retrieval problem. As we shall see it is possible that many different combinations of phase values can produce identical intensity distributions. Hence this phase retrieval problem is said to be ill-posed.
We refer back to Figure 1, where we have indicated that the contributing point sources are Fourier-transformed by an optical Fourier transform system. This type of optical system is well known; see, for example, [4]. The distribution in the Fourier plane is scaled by both the wavelength, , of the light source and the focal length of the lens. It is convenient for our purposes here to set so that we can analyze an unscaled Fourier distribution. Here this distribution is incident on an ideal recording material, allowing us to record the continuous intensity distribution. Then we perform a second unscaled inverse Fourier transform on this intensity distribution numerically. We note however that although we perform two successive Fourier operations (a forward and inverse transform), we do not arrive back at the input plane distribution since we perform the second Fourier transform operation on the intensity of the Fourier distribution not on its complex amplitude, as in a 4-f imaging system [50].
We thus begin by noting that the Fourier transform of a point source is given bywhere is the spatial location of point source “” while the parameters and refer to its phase and amplitude respectively. Hence for a sum of distinct point sources we can writeWe are however interested in the intensity distribution of the field in the Fourier plane. We remember that the intensity of a complex number is given by that number times its complex conjugate. With this relationship in mind we now write the conjugate expression for (2)where “” refers to a complex conjugate operation. We can now write an analytical expression for the intensity distribution in the Fourier plane which when simplified can be written in the following formFrom (4) we can see that we have a sum of cosines. For each given value of there is a specific spatial frequency component , and there are terms. It is easier to analyze this double summation by initially concentrating on a single spatial frequency component, i.e., for a specific value of , and we refer to such a component as . Noting the following trigonometric relationship, , we can express aswhereWe now consider the inverse Fourier transform of , which we specifically define aswhich is now in a spatial domain once again, specifically the “analysis” plane in Figure 1. Noting thatthenWe see from (9) that the complex numbers (arising from the summation of sines and cosines) multiplying the and components are complex conjugates of each other. Hence we see thatand thatWe also remember that the iterator above spans the range . Examining (10) and (11) we recognize that the only unknowns are the values for the phase parameters, . Both the initial amplitudes and the LHS of each equation are given. The LHS values are found by calculating the Fourier transform of the intensity distribution that we measure in the Fourier plane, see the “analysis plane” in Figure 1, and determining the value of the resulting complex function at the locations .
We also note that the larger the value of , the lower the number of terms from (10) and (11) that contribute to both and . When we set about trying to find from (10) and (11), we start by first setting , in which case the equations reduce to the following:and If we arbitrarily set = 0, then we can calculate from the followingwhich means that we now have values for both and .
Now setting , (10) and (11) will reduce to the following:and It is instructive to consider the next step in the procedure before we generalize the results. Setting , and substituting into (10) and (11), we arrive atand Again we see that we can calculate the LHS of (17) and (18) from an experimental measurement and by using values for and that were estimated in the previous iteration step. We summarize our theoretical results so far: (i)We set and calculate using (14)(ii)Once and are determined, we need in turn to determine and from (12) and (13)(iii)If we can find values for and , then we need to repeat the process to determine and (iv)We repeat this iterative process times to find all the unknown phase values
To proceed further we must therefore find a solution to the following coupled set of equationswhere and are the unknown phase values for a particular step in the iterative procedure. Using a software package like Mathematica it is possible to find an analytical solution for these equations; however the solution is not unique: there are 8 possible coupled values for and .
3. Algorithm Analysis
In this section we will examine some of the implications of the theoretical analysis we have just derived. First we will examine how we may reduce the number of possible solutions at each level of the iteration, reducing them from eight to two. To find a phase value for each contributing point source the algorithm must move through steps. At each step a decision on which solution to take must be made; we examine how this can be considered as a set of different possible paths through the solution space. We examine the character of the solutions produced by the algorithm. Finally we examine how to prune (reduce) the solution space by ruling out “unphysical” solutions and outline an algorithm for implementing this procedure.
3.1. From Eight Solutions to Two
We note that it is possible to rewrite the coupled equations, (19) and (20), in the following formwhere , and where we use the trigonometric relation and hence we can write the in terms of the , i.e., . Using (21) to find an expression for and subbing into (22) we can express as a quartic equation with four possible solutions [51]. For each value of there are two possible angles for that satisfy and hence eight possible solutions for and similarly eight solutions for .
This is where the difficulties with our phase retrieval approach begin in earnest. Although it is possible to find an analytical solution to the coupled equations (19) and (20), the solution has eight different possible solutions for the unknown pair: . In reality however only one of the solution pairs corresponds to the actual correct physical result. If we randomly choose from these solutions we will have a probability of 1/8 that we have chosen the physically correct answer.
Having decided on a particular choice of , we can then set and repeat the process. Hence at each step we have a 1/8 chance of guessing the physically correct solution, and for point sources we will need to make approximately guesses at each step of the algorithm. Hence the chance that we correctly guess the correct phase for each point source is , which would seem to be a vanishingly small probability. In practical experimental phase retrieval systems we may wish to find phase values where .
Fortunately, we have found from numerical simulations that although there are eight solutions possible in principal, including solutions with complex values, in practice when we substitute real physical parameters into our system of equations several of these solutions yield the same answer. We have found that at each step in the algorithm only two solutions are produced. This is the case for all simulations we have run. We do not pursue this question further here, concluding that at each step in an iterative algorithm we are able to find two solutions that are physically valid and satisfy (19) and (20).
3.2. Solution Paths
From the previous section we have found that it is possible at each stage in the iterative procedure to reduce the number of solutions to two possible answers. At each stage in the algorithm then a “binary” choice has to be made about which solution to choose. If we have contributing point sources then we must move through stages. In Figure 2, we depict how to track our progress through such a decision tree, in this example making the following choices at each stage of the algorithm, ; see the red dashed boxes in Figure 2. In general then there are possible paths in total and each possible path is identified as a unique binary number. For example, if , then there are or 8 possible solutions and if , there are or 128 solutions while if the number of possible solutions jumps to 32,768.
3.2.1. Results of Numerical Simulations
We have run numerous numerical simulations of the equations derived in Section 2 with some interesting results that we now summarize.(i)The correct phase solution is always found among the other different solutions(ii)The overwhelming majority of solutions have complex numbers for the phase and hence can be excluded as real physical solutions to the phase retrieval problem(iii)There are often several other real valued solutions in addition to the actual real solution(iv)We have run several numerical simulations with . In one run we found that there are 32,768 solutions of which are complex and the 2 remaining solutions are real valued, one of which is the correct phase solution. On a second run we found that only solutions were complex while the 6 remaining were real valued, one of which was the correct solution.
We conclude that if we had unlimited computing resources and time it would in principle be possible to calculate all possible solutions even when is very large, for example, . In practice however this calculation would result in a solution space of , which is not feasible with standard computing resources. When is relatively small and all the possible solutions can be calculated we have found that the vast majority are complex valued, and only a small number of solutions are real valued.
3.3. Pruning the Solution Space
From the preliminary analysis we have reported in the previous section, we find that the overwhelming majority of the solutions yield phase values that are complex and hence “unphysical” in nature and can be discarded. This opens up the opportunity of pruning the solution space. We will now outline a strategy that can be employed to greatly reduce the search space. This solution depends on unique binary number that identifies each particular path through the search space and we note that it will be a -bit number. As we move through the solution space we will need a variable, which we refer to as to keep track of where we are in the search tree. We will also need an iterator variable to refer to the “stage” of the algorithm; see Figure 2.
From Figure 2 we see that at the first stage, i.e., when , we identify values for and and when , we identify possible candidates for and , which we denote by and . We initially set to be all zeros and of length . For example if then initially would be set to . We also need to address the individual bits in specifically. We do this using the following notation , where refers to individual bit. We illustrate this with the following example: Consider the case when and then , , and .
The algorithm then runs as follows:(i)For a given , set and , where the length of is ,(ii)WHILE ,(1)Using Eq. (20) and Eq. (21), find and (2)IF , then = , ELSE = (3)IF are real valued then ,(4)ELSE: IF then set , ELSE , and (iii)Return
This algorithm will run through the search space checking the solutions. If a complex solution is found it will choose the alternate solution. For a given stage in the algorithm if both and are complex the algorithm will go back up the search tree, updating and choosing a different path through the search space. This algorithm will return only the first real valued solution it finds but can be modified to return all solutions that are physically significant.
When employing the “pruning” approach to search the solutions space we must still search a very large space; however the depth of the search per path will be greatly reduced; see, for example, [52]. This will significantly reduce the number of numerical operations performed. Whether the phase retrieval problem can be solved depends on the rate of branching (or the depth of the search) before a given path can be excluded from further search. From a preliminary analysis based on simulations results we find that the depth when one must search a path before an incorrect complex value is returned as the phase estimate can be modeled using statistical methods. We have found that the average path “length” approximately follows a negative exponential distribution similar to speckle intensity statistics [37, 53].
4. Conclusion
In this manuscript we have examined the problem of measuring the phase of an optical wavefield. In the Introduction, we briefly reviewed several approaches for making this measurement before concentrating on iterative phase retrieval. In these phase retrieval problems an attempt is made to estimate the phase of wavefield using usually two or more intensity distributions captured in different Linear Canonical Transform planes. An iterative FFT based algorithm is employed that tunes/modifies the phase values iteratively so that the difference between a numerically calculated intensity distribution and an experimentally captured intensity distribution is minimized. Fienup [41] compares this algorithmic approach to a gradient descent algorithm. As is often the case with such optimization algorithms, they can stagnate at local minima, producing phase estimates that are quite far from the correct physical solution. This is particularly the case when the initial guess at the phase distribution is quite far from this correct physical solution [46]. This is an important issue and different approaches for improving these algorithms are an active research topic. There is however understandably an emphasis on the numerical aspects of these algorithms.
Here we have deliberately chosen to examine a simplified version of the phase retrieval problem, assuming that we can make perfect measurements of the intensity distribution in an image and its corresponding Fourier plane. We assume that the distribution in the image plane consists of point sources with unknown phase values. Importantly we assume that we can measure the intensity distribution in the Fourier plane over an infinite extent and with an infinitely fine resolution; in short we assume that we can perfectly measure the continuous intensity distribution. We have found the following under these idealized conditions:(1)It is possible in principal to find solutions to the phase retrieval problem that satisfy the physical constraints of the problem. The algorithm for finding these solutions works iteratively but deterministically and will find among the many solutions returned the physically correct solution(2)The number of solutions is enormous; if there are contributing point sources there are possible solutions; however the overwhelming majority of these solutions are complex valued and hence “unphysical” and can be excluded as potential solution candidates(3)We can modify the algorithm so that “unphysical” solutions are excluded which we refer to as pruning. Using a pruning approach we can significantly reduce the space of possible solutions. The performance of the pruning approach can be approximately estimated using statistical techniques. A critical factor determining whether the algorithm can be computed in a realistic time is the “branching” rate or depth of search before a given solution can be excluded
It is important to also highlight the significant practical shortcomings of this theoretical analysis.
4.1. Physical Model
We model the field whose phase values we wish to recover as a set of point sources (with known amplitude) in an image plane which are separated uniformly from each other by a fixed distance, . Of course an optical wavefield is continuous in nature and so this representation is only an approximation of the real physical field. Nevertheless in standard iterative phase retrieval problems, the intensity value measured at each pixel is considered to be a point source; hence the description of the field that we use is in keeping with the standard approach in the literature.
There is however a difference in the modeling in this paper for the intensity distribution recorded in the Fourier plane. In this paper we effectively assume that we can measure the continuous Fourier plane intensity; i.e., it is as if we can record and measure the intensity over an infinite plane in the Fourier domain and with an infinitely fine spatial resolution. Making these assumptions means that we can develop convenient analytical equations relating the unknown phase values directly to experimental measurements. These analytical equations are used then to develop the iterative algorithm so that the phase values can be estimated. This contrasts with real experimental systems in that the intensity in the Fourier plane is measured only at a finite number of locations, i.e., at each pixel, and is therefore sampled. And secondly the Fourier intensity distribution is only measured over the finite extent of the CCD/CMOS array. This summarizes the theoretical differences between the analysis presented in this manuscript and the standard theoretical treatment of the iterative phase retrieval problem.
We note however that it is possible to replicate experimentally the theoretical conditions we have assumed in the manuscript. We can achieve this by making several individual measurements of the Fourier plane intensity distribution where we move the CCD/CMOS array between measurements. Then using standard digital signal processing operations we synthetically increase both the spatial resolution and the spatial extent of the measurement by stitching together the single intensity measurements.
4.2. Effects of Noise
A significant shortcoming of this theoretical approach is the effect that even small amounts of measurement noise would have on the algorithm. In every experiment there will be measurement error, for example, electronic noise in the readout of a CCD or CMOS device. If a measurement is made on a variable that has a low signal level, the relative error tends to be larger [36]. In the situation described here, we find that intensity measurements associated with higher spatial frequencies have a lower number of contributing point sources. For example, the highest spatial frequency component of the measurement has contributions from only two point sources. In general (but not always) these measurements will contain the most noise. These measurements however are used to directly estimate and , which in turn feed iteratively back into the algorithm and are repeatedly used to estimate the remaining phase values.
So we can see that to turn this theoretical approach into a practical solution would require careful experimental measurements and the development of additional signal processing algorithms to overcome unavoidable noise in these measurements. Nevertheless we hope that by approaching the phase retrieval problem from a different direction we stimulate some new interesting ideas and insights and also shed some light on the character of the solution space for the phase distributions that are recovered from intensity measurements.
4.3. Extension to 2D
In this manuscript we deliberately focused on a simple 1D theoretical model (we only consider the and spatial coordinates) so that we could examine different aspects of the phase retrieval problem. We choose to set to zero and then iteratively solve for the other phase values. Extension to 2D, where we consider both and as well as the spatial coordinates, is not necessarily straightforward and will change depending on the geometry of the contributing point sources. For example, if the contributing point sources are arranged in a square lattice, then the highest spatial frequency components arise from the two point sources that are spatially separated from each other by the largest amount, i.e., the vertices of the square. The square has four vertices; we can assign a phase value of zero to only one of those vertices, which acts as an effective reference value for the others. Sequential illumination with lines of 1D light could also be considered whereupon the 1D analysis considered here could be more easily extended.
Data Availability
This is a theoretical paper: all signal processing data can be generated from equations therein.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The author gratefully acknowledges the EU commission (H2020-TWINN-2015 HOLO, project Nr. 687328), for supporting his attendance at Summer School “Digital Holography and Related Techniques” organized in Stuttgart, 21-23 June 2017. Part of the work described in this paper results from discussions that the author had during this meeting.