Lensless Reflection Imaging of Obliquely Illuminated Objects I: Choosing a Domain for Phase Retrieval and Ptychography

Abstract

Ptychography is a lensless imaging technology that is validated from hard X-rays to terahertz spectral range. It is most attractive for extreme ultraviolet (EUV) and X-rays as optical elements are expensive and often not available. Typically, the set up involves coherently illuminated object that directs the scattered radiation normally to detector which is parallel to the object plane. Computer processing of diffraction patterns obtained when scanning the object gives the image, more precisely, the distribution of intensity and phase on its surface. However, this scheme is inefficient for EUV and X-rays due to poor reflectivity and low penetration in all materials. Reflection mode ptychography solves the problem if illumination angles do not exceed the critical angle of object material. Changing the geometry of experiment changes physical and mathematical model of image formation. Including: diffraction integral describing beam propagation from object to detector, inverse problem, optimization of object illumination angle, position and orientation of detector, choosing size and grid of coordinate and frequency computer domains. This paper considers the wavefield scattered to detector by obliquely illuminated object and determines a domain for processing of obtained scans. Solution of inverse problem with phase retrieval and resulting numerical images will be presented in the next paper.

1. Introduction

The retrieval of the field phase on the object under study, is a basis of modern lensless imaging. Ideally, the task is formulated as follows: to restore the phase of the coherent wave field from the measurements of its modulus on two parallel planes determined by the object and the detector. For this purpose, the iterative algorithm is used with calculation of the propagation of coherent radiation from the object plane to the plane of the detector. The result is the distribution of the complex field on both planes. This problem was solved for the first time, in [1] using one of the best computers at that time with ~1 megabyte memory and ~1 megaflop performance, which made it possible to analyze a large number of 32 × 32 images. The development of the ideas of this work led to the emergence of lensless imaging and ptychography as they are currently used in the visible, EUV and X-ray ranges [2,3,4,5,6,7,8,9]. Since 1970s the original algorithm [1] has been improved to be more efficient, following the progress in computer technology. New algorithms were developed to avoid the need for measurement of the intensity distribution on the object plane and to replace it by a priori information. As a result, lensless imaging methods had become applicable in a number of physics research and commercial projects [10,11,12,13]. One of the main advantages of lensless optical systems is the fundamental possibility to get rid of a problem associated with aberrations of imaging optical elements in approaching the diffraction resolution limit in microscopy, atmospheric and astronomical optics (Obviously, in order to fully realize the indicated advantage of lensless systems, it is necessary to go beyond the paraxial methods of modeling the propagation of light beams from an object to a detector).

Over the past 10–15 years, ptychography has become one of the most popular and rapidly developing lensless imaging technologies. Its scheme is shown in Figure 1. In the simplest case, ptychography consists of a registration of a sequence of interference patterns collected upon illuminating the object at partially overlapping positions produced by an object moving across a coherent beam spot, with the radiation source and detector position being fixed. Thus, the image is determined to be a distribution of the complex field amplitude (i.e., modulus and phase) on the surface of an object. It can be obtained as a result of computer processing of a set of the patterns. Using a sequence of overlapping scans instead of a single measured pattern (as in common phase retrieval method) allows one to omit a priori suggestions concerning the object optical properties.

The principles of ptychography and the very term were proposed in the works of Hoppe et al. [14,15]. However, for a long time this approach remained unnoticed due to the lack of an effective algorithm for its implementation. Such algorithm was proposed only in 2004 [16,17]. This algorithm is now known as the Ptychography Iterative Engine (PIE). However, the algorithm has a drawback—it assumes that the illuminating field is known. After 4 years, in 2008, this problem was fixed in another iterative algorithm based on the difference map approach [18,19]. In 2009 an improved algorithm PIE—Extended Ptychography Iterative Engine (ePIE)—was proposed [20].

Today ptychography is a widely used microscopy method that provides a complete (amplitude and phase) description of flat objects. Its field of application extends from terahertz [21] to hard X-rays [9]. Since there are no optical elements used, the spatial resolution is limited only by the wavelength and the detector parameters. These advantages make ptychography an attractive technique and stimulate new applications and, imaging technology (see books [22,23], review [24], and also [25]). Let us now turn to the content of this work.

The basic idea underlying ptychography and other phase reconstruction methods is the unambiguity of solving the inverse diffraction problem for the scalar wave equation; that is, unambiguous reconstruction of the complex field near the object from the known complex field near the detector. In the case of a normal location of the object, when the planes of the object and the detector are parallel and at the same time perpendicular to the beam, unambiguity is usually beyond doubt and corresponding wave field distributions are linked by the Fresnel transformations. However, in the case of an oblique position of the object to the ray (see Figure 1), the task gets more complicated. As shown in [26], the inverse problem in paraxial approximation is reduced to solving the singular integral equation.

The second idea underlying the phase reconstruction methods is the numerical calculation of the complex field near the detector from the known complex field near the object and vice versa. All currently known effective methods of such calculations are based on the fast Fourier transform [27,28,29]. For such calculations, complex fields are specified on finite lattices (digital domains). If one of the two domains is specified on either the object or the detector, then the other domain, obviously, must be dependent on it. Unfortunately, in the literature, these formulas, as well as the principles of domain selection, were not presented for an oblique object.

In the present work, within the framework of the Helmholtz equation, we propose a criterion for the uniqueness of the inverse problem solution. As well we present the formulas for calculating diffraction by a tilted object for the case of large apertures and large distances in high resolution ptychography. These formulas include the calculation of the detector domain from the object domain parameters as a development of the formulas obtained for a vertical object in our previous work [30].

2. Materials and Methods

2.1. Phase and Modulus Retrieval of Normally Located Object

Most of the modern methods of phase retrieval are, a further development of the algorithms created 30–40 years ago (see the literature cited in [20,24] and [31,32,33]). These algorithms are presented in many articles and monographs and continue to be improved. They assume that the registration of the field intensity distribution on the detector explicitly determines the modulus of the Fourier image of the field on the object in a certain range of the spatial frequencies.

The algorithms contain 4 basic operations on the complex object function $f$: (a) the Fourier transform $F$, (b) replacing its modulus with the experimentally obtained value $A$, (c) the inverse Fourier transform $F^{-1}$ and (d) the $P$ transform of the function $f$, that is a priori information about the object. The result of these transformations is a new function

$$f^{\prime }=P\left[F^{-1}\left[A\frac{F\left[f\right]}{\left|F\left[f\right]\right|}\right]\right]$$

(1)

defined in the same domain as $f$. If $f^{\prime }$ and $f$ coincide with the specified precision, then the goal is achieved. A priori information ($P$) can be the measurement of the field modulus directly behind the object or in any other plane, the shape of the object, that is, equality to zero of the field on the object plane outside a definite area, non-negative values of the field, etc.

The diagram of a typical diffraction-resolved lensless imaging microscope contains a coherent radiation source illuminating a transparent object located in the $S$ plane (see Figure 2). Unlike conventional microscopes, there are no optical elements between the object and the detector. The laser illuminating the object on the left and the computer processing the data from the detector are not shown in the diagram.

Let us show the role of the Fourier transform using the example of an angular spectrum that describes the propagation of coherent radiation from an object $S$ to a detector $S^{\prime }$ (see Figure 2) and is an exact, and consequently, nonparaxial solution of the wave equation. It defines the wave field in a half-space through spatial harmonics in the object plane:

$$\psi \left(\mathit{r}\right)={{\displaystyle \iint }}_{-\infty }^{\infty }{\phi }_0\left(\mathit{p}\right)d^2\mathit{p}exp\left\{i\mathit{p}\mathit{\rho }+i\sqrt{k^2-p^2}z\right\},\mathit{r}=\left(x,y,z\right)=\left(\mathit{\rho },z\right),$$

(2)

where $k=\frac{2\pi }{\lambda }$ is the wave number, ${\phi }_0\left(\mathit{p}\right)$ is the Fourier transform of the field distribution $\psi \left(\mathit{\rho }, z=0\right)$ in the object plane. Papers [29,34,35,36] are devoted to the analysis of this expression and its connection with other forms of the diffraction integral. Here we will consider only issues related to modeling the propagation of non-paraxial beams in the problems of image phase reconstruction and ptychography.

Assume that the detector recording the intensity is in the far zone of the object, i.e., in the Fourier plane, where:

$$z\gg \frac{a^2}{\lambda }, \rho =z\tan {\theta }_a,$$

(3)

$a$ is the size of the object, and the aperture angle ${\theta }_a$ is assumed to be fixed. Considering (2) under condition (3) and using the stationary phase method [37], it is easy to make sure that the field at the detector takes the form:

$$\psi \left(\mathit{\rho },z\right)=\frac{2\pi k}{iz}\frac{z^2}{z^2+{\rho }^2}exp\left\{ik\sqrt{z^2+{\rho }^2}\right\}{\phi }_0\left(\mathit{p}=k\frac{\rho }{\sqrt{z^2+{\mathit{\rho }}^2}}\right).$$

(4)

Expression (4), obtained from the exact solution (2) of the Helmholtz wave equation, shows that in the far zone the spatial distribution of the field $\psi \left(\mathit{\rho },z\right)$ is expressed through the Fourier transform ${\phi }_0\left(\mathit{p}\right)$ of the field on the object plane $\psi \left(\mathit{\rho },0\right)$. Moreover, when the coordinate $\rho$ on the detector changes from 0 to $\infty$, the corresponding spatial harmonic $p$ changes from 0 to $k$. In other words, there are no harmonics $p>k$ in the far field. In the general case, this expresses the fact that the limiting spatial resolution in the far zone is determined by the wavelength $\lambda =\frac{2\pi }{k}$, regardless of the type of optical system.

In practice, in phase retrieval algorithms (see review [31]), to avoid the computational complexity of exact Expression (4), a simpler one is used—the Fresnel integral, which follows from exact Expression (2) at small aperture angles ${\theta }_a$:

$$\psi \left(\mathit{\rho },z\right)=\frac{k{e}^{ikz}}{2\pi iz}{{\displaystyle \iint }}_{-\infty }^{\infty }\psi \left({\mathit{\rho }}^{\prime },0\right)d^2{\mathit{\rho }}^{\prime }exp\left\{ik\frac{{\left(\mathit{\rho }-{\mathit{\rho }}^{\prime }\right)}^2}{2z}\right\},$$

(5)

that in the far field (3) has the form:

$$\psi \left(\mathit{\rho },z\right)=\frac{2\pi k}{iz}exp\left\{ikz\left(1+\frac{{\rho }^2}{2z^2}\right)\right\}{\phi }_0\left(\mathit{p}=k\frac{\mathit{\rho }}{z}\right).$$

(6)

Let us compare the obtained expression (6) with the exact value of the field in the far zone (4). Comparison of the arguments of the Fourier transform of the initial distribution ${\phi }_0\left(\mathit{p}\right)$ shows that the use of the Fresnel integral in the far zone (6) for calculating the field at the detector is possible only for $\rho \le z$, i.e., ${\theta }_a\le 45°$ and aperture $NA=\sin {\theta }_a\le 0.71$. Otherwise, we would be talking about the harmonics of the initial distribution $p>k$ which, according to (4), are not contained in the far zone at all. Adding to this the requirement that the phase factors in (4) and (6) coincide, we obtain the condition for the applicability of the Fresnel integral in the far field (3), (see [38]):

$$\tan {\theta }_a\le 1, \frac{a^2}{\lambda }\ll z\le \frac{4}{\pi }\frac{\lambda }{{\tan }^4{\theta }_a}$$

(7)

A more detailed theoretical discussion of the issue is available in review [39].

Thus, the angular spectrum approximation (4), as well as the paraxial approximation (6), in the far zone lead to the aforementioned connection between the field intensity distribution at the detector and the Fourier transform module of the field in the object plane. Namely, the intensity on the detector is algebraically expressed through the Fourier transform of the field amplitude on the object.

2.2. Angular Spectrum and Oblique Object

Consider the propagation of reflected directional monochromatic radiation from a tilted opaque object as shown in Figure 3. The surface of the object is located on plane $S$ in the half-plane $s>0$, the detector is on the $S^{\prime }$ plane at $z=0$. The plane $S$ is inclined to the plane $x=0$ at an angle of $\theta$, and to the plane of the detector at an angle of $\pi /2-\theta$. Without loss of generality, one can suggest $0<\theta <\pi /2$. This can be achieved by choosing the appropriate coordinate system. Then, we will assume that the radiation is described by the scalar function $\psi \left(x,y,z\right)$, which is a solution of the Helmholtz equation:

$$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}+\frac{{\partial }^2\psi }{\partial z^2}+k^2\psi =0,\hspace{1em}k=\frac{2\pi }{\lambda }.$$

(8)

We also assume that the radiation propagates in the direction of increasing $z$ from left to right being a superposition

$$\psi \left(x,y,z\right)=\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{p}_x\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{p}_y\Psi \left(p_x,p_y\right){\psi }_{{\mathit{p}}_{\perp }}\left(x,y,z\right)$$

(9)

of the plane waves (“angular spectrum method”) with transverse momentum ${\mathit{p}}_{\perp }={\left(p_x,p_y\right)}^{\mathit{T}}$

$${\psi }_{{\mathit{p}}_{\perp }}\left(x,y,z\right)=e^{i{p}_{x}x+i{p}_{y}y+i\sqrt{k^2-{p_x}^2-{p_y}^2}z}.$$

(10)

One can see that the function $\Psi \left(p_x,p_y\right)$ is the Fourier transform of the function $\psi \left(x,y,0\right)$. The fact that the plane wave (10) propagates in the direction of increasing $z$ is expressed by the plus sign in front of the square root in the exponent. Let us denote unit vectors along the $x$, $y$, $z$, $s$ axes as ${\mathit{e}}_x$, ${\mathit{e}}_y$, ${\mathit{e}}_z$, ${\mathit{e}}_s$, respectively. Then the 3D wave vector is $\mathit{p}=p_x{\mathit{e}}_x+p_y{\mathit{e}}_y+\sqrt{k^2-{p_x}^2-{p_y}^2}{\mathit{e}}_z$. The plane of the object divides the space into a conditionally upper half-space, where $x\ge -z \tan \theta$ and the lower one, where $x<-z \tan \theta$. We denote the quantities referring to the upper half-space by the “$+$” sign at the bottom, and the lower half-space by the “$-$” sign at the bottom. Occurrence of the wave ${\psi }_{{\mathit{p}}_{\perp }}$ in the positive half-space is formulated as the condition:

$$Re\left(\left(\mathit{p}-\left(\mathit{p}·{\mathit{e}}_s\right){\mathit{e}}_s-\left(\mathit{p}·{\mathit{e}}_y\right){\mathit{e}}_y\right)\cdot {\mathit{e}}_z\right)\ge 0,$$

(11)

where $Re$ is the real part, since $\mathit{p}$ is complex in general case. Considering that ${\mathit{e}}_s=-\cos \theta {\mathit{e}}_z+\sin \theta {\mathit{e}}_x$ we express (11) through the angles:

$$Re\left(\sin \theta ({\tilde{p}}_x\cos \theta +\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}\sin \theta )\right)\ge 0.$$

(12)

where, for convenience, we introduced the dimensionless parameter $\tilde{\mathit{p}}=\mathit{p}/k={\left({\tilde{p}}_x,{\tilde{p}}_y,{\tilde{p}}_z\right)}^T$. Then, the symbol “~” at the top will denote quantities divided by $k$. We introduce new notation: complex angles $\vartheta$ and $\phi$ that satisfy the conditions:

$$\sin \vartheta ={\tilde{p}}_x\cos \theta +\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}\sin \theta ,$$

(13)

$$\cos \phi ={\tilde{p}}_x\sin \theta -\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}\cos \theta .$$

(14)

From Equations (13) and (14) we can express ${\tilde{p}}_x$ through the angles $\vartheta$ and $\phi$:

$${\tilde{p}}_x=\sin \vartheta \cos \theta +\cos \phi \sin \theta ,$$

(15)

Using Equation (13), the condition (12) can be written as:

$$Re\left(\sin \theta \sin \vartheta \right)\ge 0.$$

(16)

For $\left|{\mathit{p}}_{\perp }\right|\le \mathit{k}$ the vector $\tilde{\mathit{p}}$ is real. Then, the angle $\vartheta$ is equal to the angle of inclination $\tilde{\mathit{p}}$ to the plane of the object, and the angle $\phi$ is equal to the angle of inclination $\tilde{\mathit{p}}$ to the $s$ axis (see Figure 3). For all angles $0<\theta <\pi /2$, the condition (16) is equivalent to the condition $\vartheta \ge 0$, that means a wave, either leaving the object in the upper half-space for $\vartheta >0$, or surface wave moving along the object from top to bottom at $\vartheta =0$. The angular spectrum $\psi$ can be represented as an integral of plane waves (9) and as a product of the slow and fast parts:

$$\psi \left(x,y,z\right)=u\left(x,y,z\right)e^{ikz}.$$

(17)

In what follows, we will be interested in the evolution of the slow part $u\left(x,y,z\right)$ more suitable for numerical calculations. The slow part of the plane wave looks like this:

$$u_{\tilde{\mathit{p}}}\left(x,y,z\right)=e^{ik\left({\tilde{p}}_{x}x+{\tilde{p}}_{y}y+\left(\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}-1\right)z\right)}.$$

(18)

Suppose that the field of the exit wave of the object is known — it is the function $u\left(s,y\right)=u\left(x=s\sin \theta ,y,z=-s\cos \theta \right)$ defined on the plane $S$. To obtain the field on the detector $u\left(x,y,0\right)$ we will find an arbitrary harmonic $\tilde{\mathit{q}}={\left({\tilde{q}}_s,{\tilde{q}}_y\right)}^{\mathit{T}}$:

$$u_{\tilde{\mathit{q}}}\left(s,y\right)=e^{ik\left({\tilde{q}}_{s}s+{\tilde{q}}_{y}y\right)}.$$

(19)

Substituting $x=s\sin \theta ,z=-s\cos \theta$ into (18) we get:

$$u_{\tilde{\mathit{q}}}\left(s,y\right)=e^{ik\left(\left({\tilde{p}}_x\sin \theta -\left(\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}-1\right)\cos \theta \right)s+{\tilde{p}}_{y}y\right)}.$$

(20)

From Equations (19) and (20) we obtain the equations connecting ${\tilde{\mathit{p}}}_{\perp }$ and $\tilde{\mathit{q}}$:

$$\tilde{\mathit{q}}={\left({\tilde{p}}_x\sin \theta -\left(\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}-1\right)\cos \theta ,{\tilde{p}}_y\right)}^T,$$

(21)

$${\tilde{\mathit{p}}}_{\perp +,-}={\left(\pm \cos \theta \sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}+\sin \theta \left({\tilde{q}}_s-\cos \theta \right),{\tilde{q}}_y\right)}^T.$$

(22)

Using Equation (14), we rewrite Equations (21) and (22) in the form:

$$\tilde{\mathit{q}}={\left(\cos \phi +\cos \theta ,{\tilde{p}}_y\right)}^T,$$

(23)

$${\tilde{\mathit{p}}}_{\perp +,-}={\left(\pm \cos \theta \sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}+\sin \theta \cos \phi ,{\tilde{q}}_y\right)}^T,$$

(24)

It follows from Equations (24) and (15) that:

$$\sin {\vartheta }_{+,-}=\pm \sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}.$$

(25)

Thus, for each harmonic $\tilde{\mathit{q}}$ there are two solutions. One of them, according to Equations (25) and (16), corresponds to a plane wave in the upper half-space, where $\sin {\vartheta }_{+}\ge 0$, and the other to a plane wave in the lower half-space, where $\sin {\vartheta }_{-}\le 0$. Since we consider the object to be opaque, we can assume that the wave with ${\vartheta }_{+}$ corresponds to the reflected wave, and the wave in the lower half-space with ${\vartheta }_{-}$ does not exist. Therefore, choose ${\tilde{\mathit{p}}}_{+}$:

$${\tilde{\mathit{p}}}_{\perp }={\left(\cos \theta \sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}+\sin \theta \left({\tilde{q}}_s-\cos \theta \right),{\tilde{q}}_y\right)}^T.$$

(26)

It is easy to see that the imaginary part ${\tilde{\mathit{p}}}_{\perp }$ and $\sqrt{1-{\left|{\tilde{\mathit{p}}}_{\perp }\right|}^2}$ is non-negative for all $\theta \in \left[0,\pi /2\right]$. From Equation (21), $Im\left({\tilde{q}}_s\right)=0$ and Equation (26) we obtain:

$$Im\left(\sqrt{1-{\left|{\tilde{\mathit{p}}}_{\perp }\right|}^2}\right)=Im\left({\tilde{p}}_x\right)\tan \theta ,$$

(27)

$$Im\left({\tilde{p}}_x\right)=\{\begin{array}{cc}0,\hfill & {\left({\tilde{q}}_s-\cos \theta \right)}^2+{\tilde{q}}_y^2 \le 1\hfill \\ \cos \theta \sqrt{{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2-1},\hfill & {\left({\tilde{q}}_s-\cos \theta \right)}^2+{\tilde{q}}_y^2>1 \hfill \end{array}.$$

(28)

According to Equation (18), it follows that waves with complex $\tilde{\mathit{p}}$ would decay exponentially with increasing $x$ and $z$, so that the information about the corresponding frequencies $\mathit{q}$ does not reach the detector. Therefore, an unambiguous reconstruction of the $u\left(s,y\right)$ field from the $u\left(x,y,0\right)$ field needs $Im\left(\tilde{\mathit{p}}\right)=0$:

$${\left({\tilde{q}}_s-\cos \theta \right)}^2+{\tilde{q}}_y^2 \le 1.$$

(29)

Let us suppose that ${\delta }_s$ and ${\delta }_y$ are the size of the smallest detail of the object $u\left(s,y\right)$ along the directions $s$ and $y$, correspondingly (see Figure 4). Then the spectrum is restricted by the Nyquist theorem:

$$\begin{array}{c}-\frac{\pi }{{\delta }_s}\le q_s\le \frac{\pi }{{\delta }_s}\\ -\frac{\pi }{{\delta }_y}\le q_y\le \frac{\pi }{{\delta }_y}\end{array},$$

(30)

$$\begin{array}{c}-\frac{\lambda }{2{\delta }_s}\le {\tilde{q}}_s\le \frac{\lambda }{2{\delta }_s}\\ -\frac{\lambda }{2{\delta }_y}\le {\tilde{q}}_y\le \frac{\lambda }{2{\delta }_y}\end{array}.$$

(31)

Condition (29) is satisfied for all $\tilde{\mathit{q}}$ from Equation (31) when:

$$-\sqrt{1-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}+\frac{\lambda }{2{\delta }_s}\le \cos \theta \le \sqrt{1-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}-\frac{\lambda }{2{\delta }_s}.$$

(32)

In what follows, we will assume that the object is non-crystalline, and its surface changes on a scale larger than the wavelength, and $\lambda /{\delta }_s<1$, $\lambda /{\delta }_y<1$. Then, with taking into account $\theta \in \left[0,\pi /2\right]$ the expression (32) takes the form:

$$\cos \theta \le \sqrt{1-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}-\frac{\lambda }{2{\delta }_s},$$

(33)

This inequality can be rewritten by introducing the dimensionless parameter $\Phi$:

$$\Phi =\frac{1+\frac{\lambda }{2{\delta }_s}-\sqrt{1-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}{1-\cos \theta }\le 1.$$

(34)

Expanding by small parameters $\lambda /\left(2{\delta }_s\right)$ and $\lambda /\left(2{\delta }_y\right)$ and taking into account only the first terms, we can get:

$$\Phi \approx \frac{\lambda }{2{\delta }_s\left(1-\cos \theta \right)}.$$

(35)

And, for the angle $\theta \ll 1$ we obtain:

$$\Phi \approx \frac{\lambda \cos \theta }{{\delta }_s{\sin }^2\theta }.$$

(36)

Expression (36), up to a factor of 2, coincides with the definition of the $\Phi$-parameter from [40]. Thus, we can assume that Formula (34) generalizes Formula (36), obtained in [40] using the paraxial approximation $\lambda /{\delta }_s\ll 1$, $\lambda /{\delta }_y\ll 1$.

Now we can write out the final expression for the field $u$ in the form of angular spectrum:

$$\begin{array}{c}u\left(x,y,z\right)=k^2\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{\tilde{q}}_y\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{\tilde{q}}_{s}u\left({\tilde{q}}_s,{\tilde{q}}_y\right)e^{ik\left({\tilde{p}}_{x}x+{\tilde{p}}_{y}y+\left(\sqrt{1-{\left|{\tilde{p}}_{\perp }\right|}^2}-1\right)z\right)}\hfill \\ u\left({\tilde{q}}_s,{\tilde{q}}_y\right)=\frac{1}{{\left(2\pi \right)}^2}\underset{0}{\overset{\infty }{{\displaystyle \int }}}ds\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}dy u\left(s,y\right)e^{-ik({\tilde{q}}_{s}s+{\tilde{q}}_{y}y)} , \hfill \end{array}$$

(37)

where $\tilde{\mathit{p}}$ satisfies (26). From Equation (21) it follows that:

$$d\tilde{\mathit{q}}={\left(\left(\sin \theta +\frac{{\tilde{p}}_x}{\sqrt{1-{\left|{\tilde{\mathit{p}}}_{\perp }\right|}^2}}\cos \theta \right)d{\tilde{p}}_x,d{\tilde{p}}_y\right)}^T.$$

(38)

Taking the uniqueness condition (29), we obtain the field at the detector:

$$u\left(x,y,0\right)=k^2{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{p}}_y{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{p}}_x\left(\sin \theta +\frac{{\tilde{p}}_x}{\sqrt{1-{\left|{\tilde{\mathit{p}}}_{\perp }\right|}^2}}\cos \theta \right)u\left({\tilde{q}}_s,{\tilde{p}}_y\right)e^{ik\left({\tilde{p}}_{x}x+{\tilde{p}}_{y}y\right)}.$$

(39)

The real limits of integration in Equation (39) are limited by the spectrum $u\left(\tilde{\mathit{q}}\right)$ according to (31). Equation (39) means that the Fourier transforms of the object $u\left(\tilde{\mathit{q}}\right)$ and the detector $u\left(\tilde{\mathit{p}}\right)$ are related by:

$$u\left(\tilde{\mathit{p}}\right)=\left(\sin \theta +\frac{{\tilde{p}}_x}{\sqrt{1-{\left|{\tilde{p}}_{\perp }\right|}^2}}\cos \theta \right)u\left(\tilde{\mathit{q}}\right).$$

(40)

Please note that to obtain (40) we did not use the paraxial approximation. Therefore, the result is valid for the cases of long distances and large numerical apertures. Since the dependence $\tilde{\mathit{q}}\left(\tilde{\mathit{p}}\right)$ is nonlinear (see Equation (21)), the computation of $u\left(x,y,0\right)$ from $u\left(s,y\right)$ involves three "heavy" numerical operations–Fast Fourier Transform, interpolation, and Inverse Fast Fourier Transform.

2.3. Far-Field Solution

For the phase retrieval of the X-ray region of the spectrum, the far diffraction zone is usually used since the necessary condition for oversampling is satisfied in the far zone.

Let us find a solution for the far zone. Since the object is set to be in the region $s>0$, we assume that the edge of the object closest to the detector has a coordinate $s_0>0$. For convenience, instead of $x$ and $y$ coordinates, we will use the dimensionless parameters ${\alpha }_x=x/s_0$ and ${\alpha }_y=y/s_0$. Since the object begins at $s_0$ we introduce the Fourier transform of the object $u_1\left(\tilde{\mathit{q}}\right)$ beginning from $s_0$. According to (37):

$$u\left(\tilde{\mathit{q}}\right)=\left(\frac{1}{{\left(2\pi \right)}^2}{{\displaystyle \int }}_0^{\infty }ds{{\displaystyle \int }}_{-\infty }^{\infty }dy u\left(s_0+s,y\right)e^{-ik({\tilde{q}}_{s}s+{\tilde{q}}_{y}y)} \right)e^{-ik{\tilde{q}}_s{s}_0}=u_1\left(\tilde{\mathit{q}}\right)e^{-ik{\tilde{q}}_s{s}_0},$$

(41)

and the field on the detector in terms $u_1\left(\tilde{\mathit{q}}\right)$ according to Equations (37) and (41) is:

$$u\left(x,y,0\right)=k^2{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{q}}_y{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{q}}_s{u}_1\left(\tilde{\mathit{q}}\right)e^{ik\left({\tilde{p}}_x{\alpha }_x+{\tilde{p}}_y{\alpha }_y-{\tilde{q}}_s\right)s_0}.$$

(42)

The derivatives ${\tilde{p}}_x\left(\tilde{\mathit{q}}\right)$ up to the second order can be written as:

$$\begin{array}{c}\frac{\partial {\tilde{p}}_x}{\partial {\tilde{q}}_s}=\sin \theta -\frac{\cos \theta \left({\tilde{q}}_s-\cos \theta \right)}{\sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}}\\ \frac{\partial {\tilde{p}}_x}{\partial {\tilde{q}}_y}=-\frac{\cos \theta {\tilde{q}}_y}{\sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}}\end{array},$$

(43)

$$\begin{array}{c}\frac{{\partial }^2{\tilde{p}}_x}{{\partial }^2{\tilde{q}}_s}=-\frac{\cos \theta \left(1-{\tilde{q}}_y^2\right)}{{\left(\sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}\right)}^3}\\ \frac{{\partial }^2{\tilde{p}}_x}{{\partial }^2{\tilde{q}}_y}=-\frac{\cos \theta \left(1-{\left({\tilde{q}}_s-\cos \theta \right)}^2\right)}{{\left(\sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}\right)}^3}\\ \frac{{\partial }^2{\tilde{p}}_x}{\partial {\tilde{q}}_s\partial {\tilde{q}}_y}=-\frac{\cos \theta \left({\tilde{q}}_s-\cos \theta \right) {\tilde{q}}_y}{{\left(\sqrt{1-{\left({\tilde{q}}_s-\cos \theta \right)}^2-{\tilde{q}}_y^2}\right)}^3}\end{array}.$$

(44)

If the inequality (29) is fulfilled, then derivatives (43) and (44) are finite. With using the stationary phase approximation at $s_0\to \infty$ with fixed ${\alpha }_x$ and ${\alpha }_y$ we obtain:

$$\frac{\partial }{\partial {\tilde{q}}_s}\left({\tilde{p}}_x{\alpha }_x+{\tilde{p}}_y{\alpha }_y-{\tilde{q}}_s\right)=0.$$

(45)

From Equations (43) and (45) the stationary point is:

$$\begin{array}{c}{\tilde{q}}_s=\cos \theta +\frac{{\alpha }_x\sin \theta -1}{\sqrt{{\left({\alpha }_x-\sin \theta \right)}^2+{\alpha }_y^2+{\cos }^2\theta }}\\ {\tilde{q}}_y=\frac{{\alpha }_y}{\sqrt{{\left({\alpha }_x-\sin \theta \right)}^2+{\alpha }_y^2+{\cos }^2\theta }}\end{array}.$$

(46)

Then the second derivatives at the stationary point are:

$$\begin{array}{c}\frac{{\partial }^2{\tilde{p}}_x}{{\partial }^2{\tilde{q}}_s}=-\frac{1}{{\alpha }_x^3 {\cos }^2\theta }({({\alpha }_x-\sin \theta )}^2+{\cos }^2\theta )\sqrt{{({\alpha }_x-\sin \theta )}^2+{\alpha }_y^2+{\cos }^2\theta }\hfill \\ \frac{{\partial }^2{\tilde{p}}_x}{{\partial }^2{\tilde{q}}_y}=-\frac{1}{{\alpha }_x^3 {\cos }^2\theta }({\alpha }_x^2{\cos }^2\theta +{\alpha }_y^2)\sqrt{{({\alpha }_x-\sin \theta )}^2+{\alpha }_y^2+{\cos }^2\theta }\hfill \\ \frac{{\partial }^2{\tilde{p}}_x}{\partial {\tilde{q}}_s\partial {\tilde{q}}_y}=-\frac{1}{{\alpha }_x^3 {\cos }^2\theta }({\alpha }_x\sin \theta -1){\alpha }_y\sqrt{{({\alpha }_x-\sin \theta )}^2+{\alpha }_y^2+{\cos }^2\theta }\hfill \end{array}.$$

(47)

Using Equations (26) and (46), we can find ${\tilde{p}}_x$ at the stationary point:

$${\tilde{p}}_x=\frac{{\alpha }_x-\sin \theta }{\sqrt{{\left({\alpha }_x-\sin \theta \right)}^2+{\alpha }_y^2+{\cos }^2\theta }}.$$

(48)

Finally, for the field at the detector in the far zone, we obtain:

$$u\left(x,y,0\right)=\frac{2\pi k}{i{s}_0}\frac{{\alpha }_x\cos \theta }{{\left({\alpha }_x-\sin \theta \right)}^2+{\alpha }_y^2+{\cos }^2\theta }{u}_1\left(\tilde{\mathit{q}}\right)e^{ik{s}_0\left(\sqrt{{\left({\alpha }_x-\sin \theta \right)}^2+{\alpha }_y^2+{\cos }^2\theta }-\cos \theta \right)}.$$

(49)

Defining the beginning of the object on $X$ axis as $x_0=s_0\sin \theta$ and the distance from the object to the detector as $d=s_0\cos \theta$, we can get a familiar formula:

$$u\left(x,y,0\right)=\frac{2\pi k}{i}\frac{x\cos \theta }{{\left(x-x_0\right)}^2+y^2+d^2}{u}_1\left(\tilde{\mathit{q}}\right)e^{ik\left(\sqrt{{\left(x-x_0\right)}^2+y^2+d^2}-d\right)},$$

(50)

$$\mathrm{where} \tilde{\mathit{q}}={\left(\cos \theta +\frac{x\sin \theta -s_0}{\sqrt{{\left(x-x_0\right)}^2+y^2+d^2}},\frac{y}{\sqrt{{\left(x-x_0\right)}^2+y^2+d^2}}\right)}^T.$$

(51)

It is easy to show that the function $\tilde{\mathit{q}}\left(x,y\right)$ is not one-to-one. There are infinitely large values of $x$ and $y$, for which $\tilde{\mathit{q}}\left(x,y\right)$ has the same values. Therefore, for further analysis, it will be more convenient for us to use the aperture coordinates, which we define as follows:

$$R=\sqrt{{\left(x-x_0\right)}^2+y^2+d^2},$$

(52)

$$n_x=\left(x-x_0\right)/R,$$

(53)

$$n_y=y/R.$$

(54)

In aperture coordinates, Equations (50) and (51) take the form:

$$u\left(n_x,n_y,R\right)=\frac{2\pi k}{iR}\left(n_x\cos \theta +\sqrt{1-{n_x}^2-{n_y}^2}\sin \theta \right)u_1\left(\tilde{\mathit{q}}\right)e^{ikR\left(1-\sqrt{1-{n_x}^2-{n_y}^2}\right)},$$

(55)

$$\tilde{\mathit{q}}={\left(\left(1-\sqrt{1-{n_x}^2-{n_y}^2}\right)\cos \theta +n_x\sin \theta ,n_y\right)}^T.$$

(56)

In this form, the function $\tilde{\mathit{q}}\left(n_x,n_y\right)$ already one-to-one depends on $n_x$ and $n_y$.

The far-field criterion for a tilted object is also of interest, but it is outside the scope of this work. It is obvious that a use of the far-field criterion for a vertical object (3) is sufficient here.

2.4. Frequency and Spatial Domains

In order to apply Equations (39) and (50) in numerical calculations, it is necessary to find out the sizes of the frequency and spatial domains. The size of the frequency domain determines the integration limits in (39), and the size of the spatial domain determines the integration step. The spatial domain of the object is the area of the numerical definition of the function $u\left(s,y\right)$. We will treat it as a rectangular and uniform grid in coordinates $s,y$ with intervals (pixels) ${\delta }_s$, ${\delta }_y$ and with dimensions $L_s$, $L_y$ (see Figure 4). Thus, the object lies in the rectangle $\left[s_0,s_0+L_s\right]\times \left[-L_y/2,L_y/2\right]$, where $s_0$ is the coordinate of the object closest to the detector along the $S$ axis. According to Nyquist’s theorem, the frequency domain of the object is also a rectangular and uniform grid with pixels $2\pi /L_s$, $2\pi /L_y$ and dimensions $\left[-\pi /{\delta }_s,\pi /{\delta }_s\right]\times \left[-\pi /{\delta }_y,\pi /{\delta }_y\right]$. For dimensionless frequencies ${\tilde{q}}_s$ and ${\tilde{q}}_y$, the frequency domain has pixels $\lambda /L_s$, $\lambda /L_y$ and dimensions $\left[-\lambda /\left(2{\delta }_s\right),\lambda /\left(2{\delta }_s\right)\right]\times \left[-\lambda /(2{\delta }_y),\lambda /\left(2{\delta }_y\right)\right]$. Similarly, the detector domain will be treated as a rectangular uniform grid with pixels ${\delta }_x^{\prime }$, ${\delta }_y^{\prime }$ and dimensions $\left[x_{min},x_{min}+L_x^{\prime }\right]\times \left[-L_y^{\prime }/2,L_y^{\prime }/2\right]$, where $x_{min}$ is the minimum coordinate along the $X$ axis. We also introduce the following notation: $x_0=s_0\sin \theta$ is the projection of $s_0$ onto the $X$ axis, $d=s_0\cos \theta$ is the distance from the object to the detector. Further in Section 2.5, Section 2.6, Section 2.7 and Section 2.8 the parameters of the detector domain ${\delta }_x^{\prime }$, ${\delta }_y^{\prime }$, $x_{min}$, $L_x^{\prime }$ and $L_y^{\prime }$ are determined according to preset parameters of the object domain ${\delta }_s$, ${\delta }_y$, $s_0$, $L_s$ and $L_y$.

2.5. Frequency Domain of the Detector

The frequency domain $\tilde{\mathit{q}}$ is determined by inequalities (31) and, under the assumption of smoothness $u\left(s,y\right)$ at distances of the order of the wavelength, is a subset of the domain $\left[-1,1\right]\times \left[-1,1\right]$. Consider the function ${\tilde{p}}_x\left(\tilde{\mathit{q}}\right)$ from Equation (26). According to (43) the stationary point ${\tilde{p}}_x\left(\tilde{\mathit{q}}\right)$ is reached when $\tilde{\mathit{q}}$ = $\left[\cos \theta +\sin \theta ,0\right]$. It is easy to check that $\cos \theta +\sin \theta \ge 1$ at the considered angles $\theta \in \left[0,\pi /2\right]$. Therefore, the inner part of the $\tilde{\mathit{q}}$ domain does not contain the stationary point ${\tilde{p}}_x\left(\tilde{\mathit{q}}\right)$. In addition, the second derivative ${\partial }^2{\tilde{p}}_x$ is negative definite, since it follows from Equations (44), (25) and (16) that the eigenvalues ${\partial }^2{\tilde{p}}_x$ are negative and equal to $-\cos \theta /{\sin }^3\vartheta$ and $-\cos \theta /\sin \vartheta$. Thus, the minimum and maximum values of ${\tilde{p}}_x\left(\tilde{\mathit{q}}\right)$ are on the domain boundary $\tilde{q}$. It is also necessary to take into account that the radical part in (26) must be nonnegative. Therefore, there are several cases, depending on the area of intersection of the rectangular region (31) and the circle of unit radius centered at the point ${\left(cos\theta ,0\right)}^T$:

$max\left(\left|{\tilde{p}}_x\right|\right)$	Condition	(57)
$\begin{array}{c}\sin \theta \left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)\hfill \\ -\cos \theta \sqrt{1-{\left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}\hfill \end{array}$	${\left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)}^2+{\left(\frac{\lambda }{2{\delta }_y}\right)}^2\le 1$.
$\sin \theta \left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)$	${\left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)}^2+{\left(\frac{\lambda }{2{\delta }_y}\right)}^2>1,$ $\frac{\lambda }{2{\delta }_s}+\cos \theta \le 1$.
$\sin \theta$	$\frac{\lambda }{2{\delta }_s}+\cos \theta >1,$ $\theta >\arccos \left(\frac{1+\sqrt{1+8\left(1-\frac{\lambda }{2{\delta }_s}\right)}}{4}\right).$
$\cos \theta \sqrt{1-{\left(\frac{\lambda }{2{\delta }_s}-\cos \theta \right)}^2}+\sin \theta \left(\frac{\lambda }{2{\delta }_s}-\cos \theta \right)$	$\theta \le \arccos \left(\frac{1+\sqrt{1+8\left(1-\frac{\lambda }{2{\delta }_s}\right)}}{4}\right).$

Thus, Equation (57) exhausts all possible cases and allows you to analytically find the domain ${\tilde{p}}_x\in \left[-max\left(\left|{\tilde{p}}_x\right|\right),max\left(\left|{\tilde{p}}_x\right|\right)\right]$. The pixel size ${\delta }_x^{\prime }$ can be chosen based on the maximum modulus ${\tilde{p}}_x$:

${\delta }_x^{\prime }=\frac{\lambda }{2max\left(\left|{\tilde{p}}_x\right|\right)}$	Condition	(58)
$\frac{{\delta }_s}{\sin \theta \left(1+\frac{2{\delta }_s}{\lambda }\cos \theta \left(1-\sqrt{1-2\frac{\cos \theta }{si{n}^2\theta } \frac{\lambda }{2{\delta }_s}-{\left(\frac{\lambda }{2{\delta }_s}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}\right)\right)}$	${\left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)}^2+{\left(\frac{\lambda }{2{\delta }_y}\right)}^2\le 1$.
$\frac{{\delta }_s}{\sin \theta \left(1+\frac{2{\delta }_s}{\lambda }\cos \theta \right)}$	${\left(\frac{\lambda }{2{\delta }_s}+\cos \theta \right)}^2+{\left(\frac{\lambda }{2{\delta }_y}\right)}^2>1,$ $\frac{\lambda }{2{\delta }_s}+\cos \theta \le 1$.
$\frac{\lambda }{2\sin \theta }$	$\frac{\lambda }{2{\delta }_s}+\cos \theta >1,$ $\theta >\arccos \left(\frac{1+\sqrt{1+8\left(1-\frac{\lambda }{2{\delta }_s}\right)}}{4}\right).$
$\frac{\lambda }{2\left(\cos \theta \sqrt{1-{\left(\frac{\lambda }{2{\delta }_s}-\cos \theta \right)}^2}+\sin \theta \left(\frac{\lambda }{2{\delta }_s}-\cos \theta \right)\right)}$	$\theta \le \arccos \left(\frac{1+\sqrt{1+8\left(1-\frac{\lambda }{2{\delta }_s}\right)}}{4}\right).$

To determine ${\delta }_y^{\prime }$, it is enough to substitute $\theta =0$ into Equation (58). Considering that it makes no sense to choose $\lambda /\left(2{\delta }_y\right)>1$, since this no longer expands the frequency domain, we get ${\delta }_y^{\prime }={\delta }_y$.

2.6. Spatial Domain of the Detector

According to (39) and (41), the field at the detector may be expressed as:

$$u\left(x,y,0\right)=k^2{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{p}}_y{{\displaystyle \int }}_{-\infty }^{\infty }d{\tilde{p}}_x\left(\sin \theta +\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\cos \theta \right)u_1\left({\tilde{q}}_s,{\tilde{p}}_y\right)e^{ik\left({\tilde{p}}_{x}x+{\tilde{p}}_{y}y-{\tilde{q}}_s{s}_0\right)}.$$

(59)

Since the size of the region on the plane $z=0$ is $L_x^{\prime }\times L_y^{\prime }$, the steps of numerical integration in (59), according to Nyquist’s theorem, are $\left(2\pi /L_x^{\prime }\right)/k=\lambda /L_x^{\prime }$ $\times$ $\left(2\pi /L_y^{\prime }\right)/k=\lambda /L_y^{\prime }$. For numerical integration over ${\tilde{p}}_x$ to be successful, the phase change in the exponent in (59) should not exceed $2\pi$ at step $d{\tilde{p}}_x=\lambda /L_x^{\prime }$:

$$\left|k\left(d{\tilde{p}}_{x}x-d{\tilde{q}}_s{s}_0\right)\right|\le 2\pi \to \left|d{\tilde{p}}_{x}x-d{\tilde{q}}_s{s}_0\right|\le \lambda .$$

(60)

Expressing in this inequality $d{\tilde{q}}_s$ in terms of $d{\tilde{p}}_x$ according to (38), we obtain:

$$\left|x-\left(\sin \theta +\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\cos \theta \right)s_0\right|\left|d{\tilde{p}}_x\right|\le \lambda .$$

(61)

Substituting $d{\tilde{p}}_x=\lambda /L_x^{\prime }$ in (61) and using the notation $x_0$ and $d$, we get:

$$\left|\left(x-x_0\right)-\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}d\right|\le L_x^{\prime }.$$

(62)

Or, equivalently:

$$\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}d-L_x^{\prime }\le \left(x-x_0\right)\le \frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}d+L_x^{\prime }.$$

(63)

Expression (63) is valid for any $x$ and ${\tilde{p}}_x$, in particular for $x=x_{min}$, ${\tilde{p}}_x=max\left(\left|{\tilde{p}}_x\right|\right)$:

$${\left|\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\right|}_{max}d-L_x^{\prime }\le \left(x_{min}-x_0\right)$$

(64)

and for $x=x_{min}+L_x^{\prime }$, ${\tilde{p}}_x=-max\left(\left|{\tilde{p}}_x\right|\right)$:

$$\left(x_{min}-x_0+L_x^{\prime }\right)\le -{\left|\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\right|}_{max}d+L_x^{\prime }.$$

(65)

Adding inequalities (64) and (65) we obtain:

$$L_x^{\prime }\ge 2{\left|\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\right|}_{max}d=\frac{\lambda }{{\delta }_x^{\prime }} \frac{d}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(66)

Inequality (65) obviously implies the inequality:

$$x_{min}\le x_0-{\left|\frac{{\tilde{p}}_x}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}\right|}_{max}d=x_0-\frac{\lambda }{{\delta }_x^{\prime }} \frac{d/2}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(67)

Inequalities (60)–(67) are obtained under the assumption that $s_0$ is the point closest to the object. But absolutely similar inequalities can be obtained for any point of the object, including for the upper (far from the detector) edge of the object. The right-hand side in inequality (66) corresponds to the size of the diffraction spot from the point at position $s_0$. The total size of the domain $L_x^{\prime }$ is actually equal to the size of the diffraction spot from the entire object and, therefore, should be the sum of half of the diffraction spot from the near edge of the object at a distance $d$, the projection of the length of the object onto the detector $L_s\sin \theta$ and half of the diffraction spot from the far edge of the object at a distance $d+L_s\cos \theta$. As a result, you should accept:

$$L_x^{\prime }=L_s\sin \theta +\frac{\lambda }{{\delta }_x^{\prime }} \frac{\left(d+L_s\cos \theta /2\right)}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(68)

For the minimum coordinate $x_{min}$, similar reasoning leads to the fact that it is necessary to take the projection of the lower edge $x_0$ minus half of its diffraction spot. This corresponds exactly to the right-hand side of (67):

$$x_{min}=x_0-\frac{\lambda }{{\delta }_x^{\prime }} \frac{d/2}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(69)

Similarly, for $L_y^{\prime }$ we get:

$$\frac{{\tilde{p}}_y}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}d-L_y^{\prime }\le y\le \frac{{\tilde{p}}_y}{\sqrt{1-{\tilde{p}}_x^2-{\tilde{p}}_y^2}}d+L_y^{\prime }.$$

(70)

$$L_y^{\prime }=L_y+\frac{\lambda }{{\delta }_y} \frac{\left(d+L_s\cos \theta \right)}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}},$$

(71)

From Equations (68) and (69) it follows that the center of the detector domain turned out to be displaced along the $X$ axis relative to the center of the object domain by the value:

$$\Delta x_c=x_{min}+L_x^{\prime }/2-\left(s_0+L_s/2\right)\sin \theta =\frac{\lambda }{{\delta }_x^{\prime }} \frac{L_s\cos \theta /4}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(72)

We believe that it is more convenient for calculations to have a detector domain whose center has the same $x$ and $y$ coordinates as the center of the object. Therefore, let’s increase $L_x^{\prime }$ so that the new domain overlaps the old one and at the same time its center is on the same axis with the center of the object domain. To do this, obviously, you need to add $2\Delta x_c$ to $L_x^{\prime }$ and we get the final result:

$$L_x^{\prime }=L_s\sin \theta +\frac{\lambda }{{\delta }_x^{\prime }} \frac{\left(d+L_s\cos \theta \right)}{\sqrt{1-{\left(\frac{\lambda }{2{\delta }_x^{\prime }}\right)}^2-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2}}.$$

(73)

2.7. Frequency Domain of the Detector in the Far Field

In the far field, we will use aperture coordinates (52)–(54), which give an unambiguous relationship between the spatial frequency $\tilde{\mathit{q}}$ of the field on the object and the sine of the aperture angles $n_x$ and $n_y$. According to (55), the field $u\left(n_x,n_y,R\right)$ is equal to the product of the fast part

$$e^{ikR\left(1-\sqrt{1-{n_x}^2-{n_y}^2}\right)}$$

(74)

and the slow part

$$u_2\left(n_x,n_y,R\right)=\frac{2\pi k}{iR}\left(n_x\cos \theta +\sqrt{1-{n_x}^2-{n_y}^2}\sin \theta \right)u_1\left(\tilde{\mathit{q}}\right).$$

(75)

Typically, applications compute the slow part $u_2\left(n_x,n_y,R\right)$. Therefore, we define for it the frequency and spatial domains in coordinates $n_x$ and $n_y$.

Let’s designate pixels of this domain as $\delta n_x$ and $\delta n_y$, and sizes as $\left[N_{xmin},N_{xmax}\right]$ and $\left[N_{ymin},N_{ymax}\right]$. The frequency domain can be determined by noting that the change in the detector coordinates to $\delta n_x$ and $\delta n_y$ is associated with a change in the corresponding frequencies on the object $d{\tilde{q}}_s$ and $d{\tilde{q}}_y$. This connection is easy to find by differentiating (56):

$$\begin{array}{c}\left(\begin{array}{c}d{\tilde{q}}_s\\ d{\tilde{q}}_y\end{array}\right)=\\ \left(\begin{array}{cc}\sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_x+\left(1-{n_x}^2\right)\sin \theta \hfill & \sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_y-n_x n_y\sin \theta \hfill \\ -n_x n_y\hfill & 1-{n_y}^2\hfill \end{array}\right)\left(\begin{array}{c}\delta n_x\\ \delta n_y\end{array}\right).\end{array}$$

(76)

Let us reverse expression (76) and obtain the required relationship:

$$\begin{array}{c}\left(\begin{array}{c}\delta n_x\\ \delta n_y\end{array}\right)=\frac{1}{\cos \theta n_x\sqrt{1-{n_x}^2-{n_y}^2}+\left(1-{n_x}^2-{n_y}^2\right)\sin \theta }\\ \left(\begin{array}{cc}1-{n_y}^2& n_x n_y\\ n_x n_y\sin \theta -\sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_y& \sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_x+\left(1-{n_x}^2\right)\sin \theta \end{array}\right)\left(\begin{array}{c}d{\tilde{q}}_s\\ d{\tilde{q}}_y\end{array}\right).\end{array}$$

(77)

Changing by one pixel in the object domain $d{\tilde{q}}_s=\lambda /L_s$ must change at least one pixel in the detector domain $\delta n_x$ or $\delta n_y$. Therefore, at least one of the two inequalities must hold:

$$\left|\frac{\left(1-{n_y}^2\right)}{\cos \theta n_x\sqrt{1-{n_x}^2-{n_y}^2}+\left(1-{n_x}^2-{n_y}^2\right)\sin \theta }\frac{\lambda }{L_s}\right|\ge \delta n_x$$

(78)

or

$$\left|\frac{\left(n_x n_y\sin \theta -\sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_y\right)}{\cos \theta n_x\sqrt{1-{n_x}^2-{n_y}^2}+\left(1-{n_x}^2-{n_y}^2\right)\sin \theta }\frac{\lambda }{L_s}\right|\ge \delta n_y.$$

(79)

Minimizing the left-hand side of inequality (78) with the condition ${n_x}^2+{n_y}^2\le 1$, we obtain:

$$\delta n_x=\frac{2}{1+\sin \theta }\frac{\lambda }{L_s}.$$

(80)

Similarly, for $d{\tilde{q}}_y=\lambda /L_y$ we get:

$$\left|\frac{n_x n_y}{\cos \theta n_x\sqrt{1-{n_x}^2-{n_y}^2}+\left(1-{n_x}^2-{n_y}^2\right)\sin \theta }\frac{\lambda }{L_y}\right|\ge \delta n_x$$

(81)

or

$$\left|\frac{\sqrt{1-{n_x}^2-{n_y}^2}\cos \theta n_x+\left(1-{n_x}^2\right)\sin \theta }{\cos \theta n_x\sqrt{1-{n_x}^2-{n_y}^2}+\left(1-{n_x}^2-{n_y}^2\right)\sin \theta }\frac{\lambda }{L_y}\right|\ge \delta n_y.$$

(82)

Hence, minimizing the left-hand side of (82), we obtain:

$$\delta n_y=\frac{\lambda }{L_y}.$$

(83)

Finally, for the far field aperture pixel sizes, we get:

$$\delta n_x=\frac{2}{1+\sin \theta }\frac{\lambda }{L_s}$$

(84)

$$\delta n_y=\frac{\lambda }{L_y}.$$

(85)

A rectangular grating symmetrical about the optical axis in aperture coordinates is mapped onto the detector plane $S^{\prime }$ in the form of a pincushion grating with the maximum cell density in the center. Therefore, the pixels in the $S$ plane can be chosen equal:

$${\delta }_x^{\prime }=d·\delta n_x,$$

(86)

$${\delta }_y^{\prime }=d·\delta n_{\mathrm{y}}.$$

(87)

2.8. Spatial Domain of the Detector in the Far Field

Let’s now calculate the sizes of the detector domain $\left[N_{xmin},N_{xmax}\right]$ and $\left[N_{ymin},N_{ymax}\right]$. From (56) we find that $n_y={\tilde{q}}_y$, whence we immediately obtain:

$$\begin{array}{c}{N}_{ymin}=-\frac{\lambda }{2{\delta }_y}\\ N_{ymax}=\frac{\lambda }{2{\delta }_y}\end{array}.$$

(88)

The quantity $n_x$ satisfies the quadratic equation:

$$n_x^2-2\sin \theta \left({\tilde{q}}_s-\cos \theta \right)n_x+{\left({\tilde{q}}_s-\cos \theta \right)}^2-\left(1-n_y^2\right){\cos }^2\theta =0.$$

(89)

This equation has two roots, of which only the root

$$n_x=\sin \theta \left({\tilde{q}}_s-\cos \theta \right)+\cos \theta \sqrt{1-n_y^2-{\left({\tilde{q}}_s-\cos \theta \right)}^2}$$

(90)

satisfies the condition ${\tilde{q}}_s\left(n_x=0,n_y=0\right)=0$, which follows from (56). The function $n_x\left({\tilde{q}}_s,n_y\right)$ has a single maximum at the point ${\tilde{q}}_s=\cos \theta +\sin \theta$, $n_y=0$, which is not achieved in the domain (31), so we get:

$$\begin{array}{c}{N}_{xmin}=-\sin \theta \left(\frac{\lambda }{2{\delta }_s} +\cos \theta \right)+\cos \theta \sqrt{1-{\left(\frac{\lambda }{2{\delta }_y}\right)}^2-{\left(\frac{\lambda }{2{\delta }_s} +\cos \theta \right)}^2}\\ N_{xmax}=\sin \theta \left(\frac{\lambda }{2{\delta }_s} -\cos \theta \right)+\cos \theta \sqrt{1-{\left(\frac{\lambda }{2{\delta }_s} -\cos \theta \right)}^2}\end{array}.$$

(91)

The value $N_{xmin}$ becomes positive at $\theta =0$, which indicates that at small angles $\theta$ the diffraction pattern is displaced relative to the object and is completely in the upper half-space. The number of pixels in the tilt direction $\left(N_{xmax}-N_{xmin}\right)/\delta n_x$ also decreases, which indicates a decrease in the share of information about the object reaching the detector. If you want a domain that is symmetrical about an axis, you should choose $\left[-\left|N_{xmin}|,|N_{xmin}\right|\right]$ as closest spanning.

Corresponding domain size in the $S^{\prime }$ plane may be expressed as:

$$L_x^{\prime }= \frac{d·\left|N_{xmin}\right|}{\sqrt{1- {\left|N_{xmin}\right|}^2- {\left|N_{ymin}\right|}^2}}$$

(92)

$$L_y^{\prime }= \frac{d·\left|N_{ymin}\right|}{\sqrt{1- {\left|N_{xmin}\right|}^2- {\left|N_{ymin}\right|}^2}}$$

(93)

3. Results and Discussion

The results from Section 2 were implemented as a computer program that was used to validate formulas for domains. A numerical experiment (see Figure 5) was carried out in which ten thousand identical monochromatic point sources, aligned in one line with an interval of one wavelength, illuminated the plane at an angle of 35°, at a distance of 5 cm on the axis and with a numerical aperture of 0.165. The illuminated area on the plane served as an object with the parameters $\theta$ = 35°, $L_s$ = 3.84 cm, $L_y$ = 1.68 cm, $\lambda$ = 1060 nm. The detector plane was located at a distance of $d$ = 1.43 cm from the nearest edge of the object. Such a source was chosen in such a way that, on the one hand, it would satisfy the wave equation, and on the other hand, it would be narrowly directed. The latter is necessary in order for the radiation beam to be completely placed in a given numerical aperture and the illuminated spot to come to naught at the edges. In this case, one can directly compare the field of the source at the location of the detector (at a distance of 8 cm from its end) with the field calculated by Equations (39), (58), (71), (73). In this way, the accuracy with which the detector domain was chosen can be estimated. Of course, a small part of the radiation from the source still fell on the edges of the object, so the edges were forcibly smoothed in cosine at a length of 1000$\lambda$, i.e., by 1 mm. The object domain pixels were chosen ${\delta }_s=4\lambda , {\delta }_y=4\lambda$, which was enough for the correct sampling of the object but cut off the frequencies ${\tilde{q}}_s$ and ${\tilde{q}}_y$ more than 1/8 in absolute value, so the comparison was carried out near the center of the detector, where the influence of high frequencies is less than at the edges. The detector domain calculated by Equations (58), (71), (73) was ${\delta }_x^{\prime }=1.89\lambda$, ${\delta }_y^{\prime }=4\lambda$ $L_x^{\prime }=4.74$ cm, $L_y^{\prime }=2.88$ cm. Figure 6 shows the relative error in calculating the field $u\left(x,0,0\right)$ according to the formulas of Section 2.5 and Section 2.6 compared to the original field $u_S\left(x,0,0\right)$. As expected, the best accuracy is achieved near the axis (center of the domain), where the relative error averaged over the interval $x\in \left[-0.08 \mathrm{cm},0.08 \mathrm{cm}\right]$ is $7.53·{10}^{-5}$. To estimate how successful the choice of such a detector domain is, similar calculations were performed for domains with different pixel and sizes, but with the same aspect ratios ${\delta }_x^{\prime }/{\delta }_y^{\prime }$ and $L_x^{\prime }/L_y^{\prime }$. The result is shown in Figure 7, where our domain is marked with a red spot. It does not correspond to the minimum error, which is approximately 2 times smaller and is achieved for pixels $1.3{\delta }_x^{\prime }$, $1.3{\delta }_y^{\prime }$ and sizes $0.8L_x^{\prime }$, $0.8L_y^{\prime }$, but nevertheless close to this. Similar comparisons were made for the far zone (see Figure 8 and Figure 9). Here, the average relative error was calculated in the range $n_x\in \left[-0.1,0.1\right]$ and amounted to $5.34·{10}^{-5}$. At the same time, the minimum is 32% less and is achieved at $2.0{\delta }_x^{\prime }$, $2.0{\delta }_y^{\prime }$ and $0.6L_x^{\prime }$, $0.6L_y^{\prime }$. Thus, the domain suggested in Section 2.5, Section 2.6, Section 2.7 and Section 2.8 is more of a compromise than is absolutely the right choice. It provides reasonable accuracy while reasonably saving computer memory.

4. Conclusions

The paper presents a model of lensfree image formation in reflection schemes with oblique object illumination. The relation between the wavefields on the object and the detector is obtained with the help of angular spectrum method and scalar wave equation.

For the first time, formulas are obtained that relate the sizes and steps of the grid of the computational domain with the geometry of experiment with inclined object.

The accuracy of the proposed domain is assessed by comparing it with the exact solution. It is shown that the specified domain provides an accuracy comparable in order of magnitude with the maximum possible, without leading to an excessive use of computer memory.

In the far field, the formulas obtained can be simplified and the algorithms become faster. More details will be presented in the next paper, where inverse problems in phase retrieval and ptychography will be considered.