Application of the Five-Step Phase-Shifting Method in Reflective Ghost Imaging for Efficient Phase Reconstruction

Abstract

The conventional approach to phase reconstruction in Reflective Ghost Imaging (RGI) typically involves the introduction of three reference screens into the reference path, deeming the Fourier transform step indispensable. However, this method introduces complexity to the system and raises concerns regarding potential errors in phase retrieval. In response to these challenges, we advocate for adopting the Five-Step Phase-Shifting (FSPS) method in the RGI system. This method presents two key advantages over traditional approaches: (1) It streamlines the phase reconstruction process by eliminating the requirement for a Fourier inverse transform. (2) It avoids the need to insert objects into the reference optical path, simplifying the computation of reference optical path intensity and enabling seamless application to Computational Ghost Imaging (CGI), overcoming the constraints of Dual-Arm Ghost Imaging (DAGI). We substantiate the theoretical proposition through numerical simulations involving two intricate objects. Furthermore, our discussion delves into exploring the influence of varying reflective angles on the phase reconstruction performance.

1. Introduction

Ghost Imaging (GI) has garnered substantial attention in recent years due to its nonlocal probing and single-pixel detecting properties. It reconstructs the image of an unknown object by measuring the intensity correlation between test and reference path detectors [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. Recent years have witnessed significant advancements in signal processing through the application of deep learning techniques [19]. These methods leverage the capabilities of neural networks to comprehend intricate patterns and relationships, providing a potent tool for further enhancing the capabilities of ghost imaging across various applications [20,21]. On the one hand, GI technology excels in low-light imaging, demonstrating proficiency in producing clear images from weak light signals. This capability lends itself to a wide range of applications, including remote sensing [22,23], watermark technology [24], medical imaging [25], and cryptography [26,27]. On the other hand, GI can handle not only transmissive objects [1,16,28,29,30,31,32,33,34,35,36,37,38,39] but also reflective ones [40,41,42,43,44,45,46,47,48]. Numerous reports indicate the successful reconstruction of both amplitude and phase distribution of transmissive objects [16,28,29,30,31,32,33,34,35,36,37,38,39] in Transmissive Ghost Imaging (TGI).

Obtaining phase information about objects is crucial in various fields, especially in optics [49], electronics [50], communications [51], medicine [52], and other scientific and engineering domains. For instance, in medicine, phase information plays a vital role in various imaging techniques, including magnetic resonance imaging and optical coherence tomography. Analyzing the phase information allows doctors to gain insights into tissue structures, vascular distribution, and abnormalities for diagnosis and treatment planning.

The commonly employed GI technique is used to determine the reflectivity of a reflective object, but it leaves the phase distribution of the reflection coefficient unknown. In recent years, efforts have been made to reconstruct the phase distribution information in Reflective Ghost Imaging (RGI) [46]. In [46], Chen successfully reconstructed the reflected object’s phase distribution by inserting reference screens into the reference path and employing the Three-Step Phase-Shifting (TSPS) method. However, the inclusion of a Fourier transform step in this method is necessary for reconstructing phase information, potentially introducing a phase retrieval error. Moreover, the incorporation of a reference screen in the reference arm restricts the applicability of this approach to computational ghost imaging (CGI) [3,4,5,6,7] and confines its use to Dual-Arm Ghost Imaging (DAGI) applications.

In a recent contribution, Chen introduced an innovative approach that utilizes the Five-Step Phase-Shifting (FSPS) method for reconstructing the phase distribution of a complex transmissive object in TGI [34]. This method offers distinct advantages by eliminating the need for iterative algorithms [35,36,37], a complex GI system [30,53], entangled photon pairs [39], or Fourier transform steps [28,29,33,46,54]. Additionally, the FSPS method transcends the limitations of DAGI and can be seamlessly applied to CGI.

Building upon these advancements, we spontaneously embarked on exploring the potential application of the novel FSPS method to RGI for obtaining the phase distribution of reflective objects. While the FSPS method was initially proposed in the context of scenarios involving TGI, our detailed theoretical derivations focused on RGI scenarios, revealing its direct applicability to enhance phase reconstruction in RGI.

In this paper, we take the initial step of providing a detailed theoretical derivation of the FSPS method in RGI. It is worth noting that, although we utilize a DAGI framework for clarity, the current framework can be seamlessly adapted for CGI as long as the reference arm contains no objects. Subsequently, we investigate the impact of various incident and reflective angles on this phase reconstruction method through numerical simulations. Our deliberate selection of the most suitable incident and reflective angles aims to showcase the robustness and efficiency of the proposed phase retrieval approach. In summary, our work bridges the gap between the FSPS method initially proposed for TGI and its direct application to enhance phase reconstruction in RGI. The adaptability of the FSPS method and its robust performance in diverse scenarios highlight its potential as a valuable tool in advancing RGI techniques.

2. Model and Theory

We present our generalized Reflective Ghost Imaging (RGI) scheme in Figure 1. Notably, it closely resembles the conventional RGI system. The classical source undergoes division into two beams facilitated by the beam splitter (BS). Subsequently, these two beams traverse distinct paths: the reference path and the test path. Within the test path, an unknown reflective object is introduced, with $d_2$ denoting the distance from the object to the test detector $D_t$, and $d_1$ representing the distance from the source to the object. Simultaneously, the reference path remains unaffected by the object and incorporates a high-resolution sensor $D_r$. The distance between the source and the detector $D_r$ is designated as $d_0$.

Theoretically, the reflective object can be described by the speckle model, with its random surface height assumed to follow Gaussian statistics. As discussed in Refs. [44,55], the relationship between the reflected fields $E_o$ at the surface of the object and the incident field $E_i$ is expressed as follows:

$$\begin{array}{c}\hfill E_o\left(u\right)=E_i\left(u\right)r\left(u\right)exp\left(j\varphi \left(u\right)\right).\end{array}$$

(1)

Here, $r\left(u\right)$ represents the reflection coefficient of the object, and $\varphi \left(u\right)$ signifies the phase delay proportional to the random surface height of the object $h\left(u\right)$:

$$\begin{array}{c}\hfill \varphi \left(u\right)=k(-\overrightarrow{i}·\overrightarrow{n}+\overrightarrow{o}·\overrightarrow{n})h\left(u\right)\end{array}$$

(2)

where k represents the wave number of the incident light, and the dot product $\overrightarrow{o}·\overrightarrow{n}$ is the cosine of the angle between the unit vectors $\overrightarrow{o}$ and $\overrightarrow{n}$.

Using Equation (1), the field $E_t\left(x_t\right)$ at the test detector $D_t$ can be calculated as [56]

$$\begin{array}{ccc}\hfill E_t\left(x_t\right)& =& \frac{-1}{\lambda \sqrt{d_1{d}_2}}\int dud{x}_1{E}_s\left(x_1\right)r\left(u\right)exp\left[i\varphi \left(u\right)\right]\hfill \\ & & \times exp\left[-\frac{i\pi }{\lambda d_1}{\left(x_1-ucos{\theta }_i\right)}^2\right]\hfill \\ & & \times exp\left[-\frac{i\pi }{\lambda d_2}{\left(ucos{\theta }_o-x_t\right)}^2\right]\hfill \end{array}$$

(3)

where ${\theta }_i$ and ${\theta }_o$ represent the incident angle and reflective angle, respectively. Similarly, we obtain the field $E_r\left(x_r\right)$ as follows:

$$\begin{array}{ccc}\hfill E_r\left(x_r\right)& =& \frac{1}{\sqrt{i\lambda d_0}}\int d{x}_2{E}_s\left(x_2\right)exp\left[-\frac{i\pi }{\lambda d_0}{(x_2-x_r)}^2\right].\hfill \end{array}$$

(4)

Combined with classical optical coherent theory [56] and the notations in [44], and considering that a complex circular Gaussian random process can model the field fluctuations of a classical light source with zero mean [57], we have

$$\begin{array}{ccc}\hfill \langle E_s\left(x_1\right)E_s^{*}\left(x_1^{\prime }\right)E_s\left(x_2\right)E_s^{*}\left(x_2^{\prime }\right)\rangle & =& \langle E_s\left(x_1\right)E_s^{*}\left(x_2^{\prime }\right)\rangle \langle E_s^{*}\left(x_1^{\prime }\right)E_s\left(x_2\right)\rangle \hfill \\ & & +\langle E_s\left(x_1\right)E_s^{*}\left(x_1^{\prime }\right)\rangle \langle E_s\left(x_2\right)E_s^{*}\left(x_2^{\prime }\right)\rangle \hfill \end{array}$$

(5)

where $\langle \rangle$ denotes the ensemble average, the correlation of intensity fluctuations between the test and reference detectors can be calculated as follows:

$$\begin{array}{ccc}\hfill G(x_t,x_r)& =& \langle E_r\left(x_r\right)E_r^{*}\left(x_r\right)E_t\left(x_t\right)E_t^{*}\left(x_t\right)\rangle -\langle E_t\left(x_t\right)E_t^{*}\left(x_t\right)\rangle \langle E_r\left(x_r\right)E_r^{*}\left(x_r\right)\rangle \hfill \\ & =& \frac{1}{{\lambda }^3{d}_0{d}_1{d}_2}\int dud{u}^{{}^{\prime }}\langle r\left(u\right)r^{*}\left(u^{{}^{\prime }}\right)exp[i\varphi \left(u\right)-i\varphi \left(u^{{}^{\prime }}\right)]\rangle \hfill \\ & & \times \int d{x}_{1}d{x}_1^{{}^{\prime }}d{x}_{2}d{x}_2^{{}^{\prime }}\langle E_s^{*}\left(x_1^{{}^{\prime }}\right)E_s\left(x_2\right)\rangle \langle E_s^{*}\left(x_2^{{}^{\prime }}\right)E_s\left(x_1\right)\rangle \hfill \\ & & \times exp\left\{\frac{i\pi }{\lambda d_1}[{(x_1^{{}^{\prime }}-u^{{}^{\prime }}cos{\theta }_i)}^2-{(x_1-ucos{\theta }_i)}^2]\right\}\rangle \hfill \\ & & \times exp\left\{\frac{i\pi }{\lambda d_2}[{(u^{{}^{\prime }}cos{\theta }_o-x_t)}^2-{(ucos{\theta }_o-x_t)}^2]\right\}\hfill \\ & & \times exp\left\{\frac{i\pi }{\lambda d_0}[{(x_2^{{}^{\prime }}-x_r)}^2-{(x_2-x_r)}^2]\right\}.\hfill \end{array}$$

(6)

In TGI, only one stochastic fluctuation originates from the source field. However, in RGI, the scenario is more intricate as it involves two independent types of stochastic fluctuations: one from the source and the other from the object itself [46]. Let us define the stochastic fluctuation from the object as $R(u,u^{\prime })=\langle r\left(u\right)r^{*}\left(u^{\prime }\right)exp[i\varphi \left(u\right)-i\varphi \left(u^{\prime }\right)]\rangle$, and it is independent of the source fluctuation. Mathematically, $R(u,u^{\prime })=r\left(u\right)r^{*}\left(u^{\prime }\right)exp\{-{\sigma }_{\varphi }^2[1-exp(-{(u-u^{\prime })}^2/{\ell }_c^2)]\}$ [58], where $r\left(u\right)$ represents the reflection coefficient of the object, ${\ell }_c$ is the surface correlation length, and ${\sigma }_{\varphi }^2$ is the variance of the phase related to the variance of the surface height ${\sigma }_h^2$,

$$\begin{array}{c}\hfill {\sigma }_{\varphi }^2=\left[k(-\overrightarrow{i}·\overrightarrow{n}+\overrightarrow{o}·\overrightarrow{n})\right]{\sigma }_h^2.\end{array}$$

(7)

Suppose the source is entirely incoherent; we have

$$\begin{array}{ccc}& & \langle E_s^{*}\left(x_2^{{}^{\prime }}\right)E_s\left(x_1\right)\rangle =f\left(x_1\right)\mathrm{rect}\left(\frac{x_1}{s_e}\right)\delta (x_1-x_2^{{}^{\prime }})\hfill \\ & & \langle E_s^{*}\left(x_1^{{}^{\prime }}\right)E_s\left(x_2\right)\rangle =f\left(x_2\right)\mathrm{rect}\left(\frac{x_2}{s_e}\right)\delta (x_2-x_1^{{}^{\prime }})\hfill \end{array}$$

(8)

in which $f\left(x_1\right)$ represents the intensity distribution of the source, $\mathrm{rect}\left(x\right)$ is the rectangular function, and $s_e$ denotes the size of the incoherent source.

In implementing GI, we set $d_0=d_1$, and position a point detector in the test path. The expression for the GI formula can be articulated as follows:

$$\begin{array}{ccc}\hfill I\left(x_r\right)& =& G(x_t=0,x_r)\hfill \\ & =& \frac{1}{{\lambda }^3{d}_0^2{d}_2}\int d{x}_{1}d{x}_{2}dud{u}^{{}^{\prime }}R(u,u^{{}^{\prime }})f\left(x_1\right)f\left(x_2\right)\mathrm{rect}\left(\frac{x_1}{s_e}\right)\mathrm{rect}\left(\frac{x_2}{s_e}\right)\hfill \\ & & \times exp\left\{\frac{i\pi }{\lambda d_0}[{(x_2-u^{{}^{\prime }}cos{\theta }_i)}^2-{(x_1-ucos{\theta }_i)}^2+{(x_1-x_r)}^2-{(x_2-x_r)}^2]\right\}\hfill \\ & & \times exp\left[\frac{i\pi {cos}^2{\theta }_o}{\lambda d_2}({u^{{}^{\prime }}}^2-u^2)\right].\hfill \end{array}$$

(9)

To simplify the formula, we introduce the following symbols:

$$\begin{array}{c}\hfill k_0=\frac{i\pi }{\lambda d_0},k_2=\frac{i\pi {cos}^2{\theta }_o}{\lambda d_2{cos}^2{\theta }_i},A=\frac{1}{{\lambda }^3{d}_0^2{d}_2}.\end{array}$$

Moreover, we substitute variables with $y=ucos{\theta }_i$ and $y^{\prime }=u^{\prime }cos{\theta }_i$. Assuming the incoherent source size is sufficiently large, the formula can be presented as:

$$\begin{array}{ccc}\hfill I\left(x_r\right)& =& \frac{A}{{cos}^2{\theta }_i}\int d{x}_{1}d{x}_{2}f\left(x_1\right)f\left(x_2\right)\hfill \\ & & \times dyd{y}^{{}^{\prime }}R(\frac{y}{cos{\theta }_i},\frac{y^{{}^{\prime }}}{cos{\theta }_i})exp\left[k_2({y^{{}^{\prime }}}^2-y^2)\right]\hfill \\ & & \times exp\left\{k_0[{(x_2-y^{{}^{\prime }})}^2-{(x_1-y)}^2]\right\}\hfill \\ & & \times exp\left\{k_0[{(x_1-x_r)}^2-{(x_2-x_r)}^2]\right\}\hfill \\ & =& \frac{A}{{cos}^2{\theta }_i}\int dyd{y}^{{}^{\prime }}R(\frac{y}{cos{\theta }_i},\frac{y^{{}^{\prime }}}{cos{\theta }_i})\hfill \\ & & \times exp\left[(k_0+k_2)({y^{{}^{\prime }}}^2-y^2)\right]\hfill \\ & & \times d{x}_{1}f\left(x_1\right)exp\left[2k_0{x}_1(y-x_r)\right]\hfill \\ & & \times d{x}_{2}f\left(x_2\right)exp\left[2k_0{x}_2(x_r-y^{{}^{\prime }})\right].\hfill \end{array}$$

(10)

Especially, when ${\ell }_c$ is significantly large, we obtain $R(u,u^{\prime })=r\left(u\right)r^{*}\left(u^{\prime }\right)$, resulting in the GI pattern as:

$$\begin{array}{ccc}\hfill I\left(x_r\right)& =& \frac{A}{{cos}^2{\theta }_i}{|\int dyr\left(\frac{y}{cos{\theta }_i}\right)exp[-(k_0+k_2)y^2]d{x}_{1}f\left(x_1\right)exp\left[2k_0{x}_1(y-x_r)\right]|}^2\hfill \\ & =& \frac{A}{{cos}^2{\theta }_i}{\left|\int dyF\left(\frac{x_r-y}{i\pi /k_0}\right)r\left(\frac{y}{cos{\theta }_i}\right)exp[-(k_0+k_2)y^2]\right|}^2\hfill \\ & =& \frac{A}{{cos}^2{\theta }_i}{\left|F\left(\frac{x_r}{i\pi /k_0}\right)\otimes \left(r\left(\frac{x_r}{cos{\theta }_i}\right)exp[-(k_0+k_2)x_r^2]\right)\right|}^2\hfill \end{array}$$

(11a)

$$\begin{array}{ccc}& =& \frac{1}{{\lambda }^3{d}_0^2{d}_2{cos}^2{\theta }_i}{\left|F\left(\frac{x_r}{\lambda d_0}\right)\otimes \left(r\left(\frac{x_r}{cos{\theta }_i}\right)exp[-\frac{i\pi }{\lambda }(\frac{1}{d_0}+\frac{1}{d_2}\times \frac{{cos}^2{\theta }_o}{{cos}^2{\theta }_i})x_r^2]\right)\right|}^2\hfill \end{array}$$

(11b)

in which ⊗ denotes convolution, and $F(.)$ represents the Fourier transform of $f(.)$. It is evident that Equation (11b) bears a strong resemblance to Equation (7) introduced in [34], suggesting that the FSPS method, originally proposed for TGI in [34], could be extended to RGI.

For simplicity, we substitute the object $r_e\left(x_r\right)$ for $r\left(\frac{x_r}{cos{\theta }_i}\right)exp\left[-\frac{i\pi }{\lambda }\left(\frac{1}{d_0}+\frac{1}{d_2}\times \frac{{cos}^2{\theta }_o}{{cos}^2{\theta }_i}\right)x_r^2\right]$. Consequently, we obtain:

$$\begin{array}{c}\hfill r_e\left(x_r\right)=r\left(\frac{x_r}{cos{\theta }_i}\right)exp[-\frac{i\pi }{\lambda }(\frac{1}{d_0}+\frac{1}{d_2}\times \frac{{cos}^2{\theta }_o}{{cos}^2{\theta }_i})x_r^2].\end{array}$$

(12)

By acquiring the phase and amplitude details of the object $r_e\left(x_r\right)$, it becomes evident that we can subsequently reconstruct the phase and amplitude distributions of $r\left(x_r\right)$ using Equation (12).

Utilizing Equation (12), we can represent Equation (11b) as follows:

$$\begin{array}{c}\hfill I\left(x_r\right)=\frac{1}{{\lambda }^3{d}_0^2{d}_2{cos}^2{\theta }_i}{\left|F\left(\frac{x_r}{\lambda d_0}\right)\otimes r_e\left(x_r\right)\right|}^2.\end{array}$$

(13)

It is evident that Equation (13) is identical to Equation (9) derived in [34], reinforcing the idea that the FSPS method proposed in [34] for TGI can be directly applied to RGI. Thus, with the FSPS method [34], one can directly obtain:

$$\begin{array}{ccc}\hfill H\left(x_r\right)& =& G_1{I}_1\left(x_r\right)+G_2{I}_2\left(x_r\right)+G_3{I}_3\left(x_r\right)+G_4{I}_4\left(x_r\right)+G_5{I}_5\left(x_r\right)\hfill \end{array}$$

(14a)

$$\begin{array}{ccc}& =& r_e(x_r-\epsilon \lambda d_0)r_e^{*}(x_r+\epsilon \lambda d_0)\hfill \end{array}$$

(14b)

where $G_1=\frac{(2\sqrt{2}-2)(1+i)}{2\sqrt{2}-4}$, $G_2=\frac{[-2\sqrt{2}+i(2\sqrt{2}-4)]}{2\sqrt{2}-4}$, $G_3=\frac{4}{2\sqrt{2}-4}$, $G_4=\frac{[-2\sqrt{2}-i(2\sqrt{2}-4)]}{2\sqrt{2}-4}$, $G_5=\frac{(2\sqrt{2}-2)(1-i)}{2\sqrt{2}-4}$, and $I_m\left(x_r\right)=\frac{1}{{\lambda }^3{d}_0^2{d}_2{cos}^2{\theta }_i}{\left|r_e\left(x_r\right)+u_m{r}_e(x_r-\epsilon \lambda d_0)+u_m^{*}{r}_e(x_r+\epsilon \lambda d_0)\right|}^2$, with $m=1,2,3,4,5$ and $u_1=\frac{1}{2}$, $u_2=\frac{1+i}{2\sqrt{2}}$, $u_3=\frac{i}{2}$, $u_4=\frac{-1+i}{2\sqrt{2}}$, $u_5=-\frac{1}{2}$. Here, $\epsilon$ denotes a real constant.

Then we have

$$\begin{array}{c}\hfill {\Phi }_H\left(x_r\right)={\Phi }_{r_e}(x_r-\epsilon \lambda d_0)-{\Phi }_{r_e}(x_r+\epsilon \lambda d_0)\end{array}$$

(15)

where ${\Phi }_H\left(x_r\right)$ represents the phase of $H\left(x_r\right)$, and ${\Phi }_{r_e}\left(x_r\right)$ represents the phase of $r_e\left(x_r\right)$, the quantitative reconstruction of the phase of $r_e\left(x_r\right)$ from the phase information of $H\left(x_r\right)$ can be achieved using Equation (15). The specific steps are as follows: (1) Assume that the phase at $x_r=0$ is zero, i.e., ${\Phi }_{r_e}\left(0\right)=0$. (2) By applying Equation (15), one can determine the values of ${\Phi }_{r_e}\left(0\right)$, $\pm {\Phi }_{r_e}\left(2\epsilon \lambda d_0\right)$, $\pm {\Phi }_{r_e}\left(4\epsilon \lambda d_0\right)$, and so forth.

Upon the successful reconstruction of the phase information of $r_e\left(x_r\right)$ from the phase information of $H\left(x_r\right)$, the phase information of $r\left(\frac{x_r}{cos{\theta }_i}\right)$ can be obtained using Equation (12):

$$\begin{array}{c}\hfill {\Phi }_r\left(\frac{x_r}{cos{\theta }_i}\right)={\Phi }_{r_e}\left(x_r\right)+\frac{\pi }{\lambda }\left(\frac{1}{d_0}+\frac{1}{d_2}\times \frac{{cos}^2{\theta }_o}{{cos}^2{\theta }_i}\right)x_r^2\end{array}$$

(16)

The amplitude of $r_e\left(x_r\right)$ can be determined using the standard RGI scheme without the designed shaped source. By applying Equation (12), we can directly infer that the amplitude of $r_e\left(x_r\right)$ is equivalent to the amplitude information of $r\left(\frac{x_r}{cos{\theta }_i}\right)$. Therefore,

$$\begin{array}{c}\hfill \left|r\left(\frac{x_r}{cos{\theta }_i}\right)\right|=\left|r_e\left(x_r\right)\right|.\end{array}$$

(17)

With Equations (16) and (17), we can ascertain the amplitude and phase of $r\left(\frac{x_r}{cos{\theta }_i}\right)$. Subsequently, the phase and amplitude information of $r\left(x_r\right)$ can be directly reconstructed from $r\left(\frac{x_r}{cos{\theta }_i}\right)$ through coordinate transformation. Assuming a zero phase at $r_e\left(x_r\right)=0$ denoted as ${\varphi }_{r_e}\left(0\right)=0$, the reconstructed phase exhibits a constant value difference from the actual phase. Nevertheless, the absolute phase holds little significance, given that the relative phase distribution remains unchanged.

3. Numerical Simulations

In the subsequent discussions, we verify the effectiveness of our RGI scheme by employing two types of complex reflective objects: a Reflective Double-Slit (RDS) and a Reflective Double-Slit Gaussian Phase Plate (RDSGPP). In our simulations, we configure the transverse size of the source as $D_s=10\mathrm{mm}$, the wavelength as $\lambda =628\mathrm{nm}$, and the distances $d_1=400\mathrm{mm}$ and $d_2=200\mathrm{mm}$. The Charge-Coupled Device (CCD) resolution is ${\Delta }_{x_r}=8.3\mathsf{\mu }\mathrm{m}$, the sample number is $M=320$, and $\epsilon$ is specified as $33.041{\mathrm{m}}^{-1}$. Here, let us elaborate on why the value of $\epsilon$ is chosen as $33.041{\mathrm{m}}^{-1}$. The reason can be traced back to Equation (15). The discretized form of Equation (15) is expressed as follows:

$$\begin{array}{ccc}\hfill {\Phi }_H(M\Delta x_r)& =& {\Phi }_{r_e}(M\Delta x_r-\Delta x_r)-{\Phi }_{r_e}(M\Delta x_r+\Delta x_r)\hfill \\ & =& {\Phi }_{r_e}[(M-1)\Delta x_r]-{\Phi }_{r_e}[(M+1)\Delta x_r]\hfill \end{array}$$

(18)

where $\Delta x_r=\epsilon \lambda d_0$, and $\Delta x_r$ represents the pixel size of the CCD ($\Delta x_r=8.3\mathsf{\mu }\mathrm{m}$). Additionally, $x_r$ must be divisible by $\Delta x_r$ (i.e., $x_r/\Delta x_r=M$, where M is an integer, and here we set $M=320$). This discrete representation proves to be more convenient for subsequent signal processing in our simulated experiments. According to Equation (18), we can extract the values of ${\Phi }_{r_e}$ from ${\Phi }_H$ using a recursive algorithm. Thus, we determine the value of $\epsilon$ as:

$$\epsilon =\frac{\Delta x_r}{\lambda d_0}=\frac{8.3\mathsf{\mu }\mathrm{m}}{628\mathrm{nm}\times 400\mathrm{mm}}=33.041{\mathrm{m}}^{-1}.$$

Here, we also outline the steps of our simulation experiment as follows:

1. Perform multidimensional integration on Equation (11a) to calculate the values of $I_1\left(x_r\right)$, $I_2\left(x_r\right)$, $I_3\left(x_r\right)$, $I_4\left(x_r\right)$, $I_5\left(x_r\right)$, and $I_0\left(x_r\right)$. These values can be determined by setting the parameter u in the function $f\left(x_1\right)=1+u{e}^{i2\pi \epsilon x_1}+u^{*}{e}^{-i2\pi \epsilon x_1}$ [34] to $u=\frac{1}{2}$, $u=\frac{1+i}{2\sqrt{2}}$, $u=\frac{i}{2}$, $u=\frac{-1+i}{2\sqrt{2}}$, $u=-\frac{1}{2}$, and $u=0$. Obtaining the value of $I_0\left(x_r\right)$ is crucial for extracting the amplitude information of the object: setting u to 0 transforms it into a conventional RGI mode. Taking the square root of $I_0\left(x_r\right)$ yields the amplitude information $\left|r_e\left(x_r\right)\right|$ of $r_e\left(x_r\right)$.

2. Substitute the values of $I_1\left(x_r\right)$, $I_2\left(x_r\right)$, $I_3\left(x_r\right)$, $I_4\left(x_r\right)$, and $I_5\left(x_r\right)$ obtained in the first step into Equation (14a), resulting in $H\left(x_r\right)$. Then, obtain the phase value ${\Phi }_H\left(x_r\right)$ from $H\left(x_r\right)$. By further applying Equation (18) and using a recursive algorithm, we can extract the phase values of ${\Phi }_{r_e}\left(x_r\right)$.

3. Substitute the obtained values of ${\Phi }_{r_e}\left(x_r\right)$ from the second step into Equation (16), successfully obtaining the phase information ${\Phi }_r\left(\frac{x_r}{cos{\theta }_i}\right)$. Then, applying Equation (17) along with the amplitude $\left|r_e\left(x_r\right)\right|$ obtained in the first step yields the value of $\left|r\left(\frac{x_r}{cos{\theta }_i}\right)\right|$.

4. Transform the obtained phase values ${\Phi }_r\left(\frac{x_r}{cos{\theta }_i}\right)$ and amplitude values $\left|r\left(\frac{x_r}{cos{\theta }_i}\right)\right|$ using coordinate transformations to obtain the values of ${\Phi }_r\left(x_r\right)$ and $\left|r\left(x_r\right)\right|$. Thus, the reconstruction of the phase and amplitude of the reflective object is completed.

These steps are further presented in the form of a flowchart in Figure 2.

To investigate the influence of different incident angles ${\theta }_i$ and reflective angles ${\theta }_o$ on our phase retrieval method, we exemplify using the RDSGPP and simulate the reconstructed results for nine distinct angle combinations presented in Figure 3. The reflection coefficient of the RDSGPP is defined as:

$$r_1\left(u\right)=\left\{\begin{array}{cc}{\rho }_1{e}^{-\frac{i{u}^2}{{{\tau }_1}^2}},\hfill & \mathrm{if}-\frac{2w+d}{2}\le u\le -\frac{d}{2},\hfill \\ {\rho }_2{e}^{-\frac{i{u}^2}{{{\tau }_2}^2}},\hfill & \mathrm{if}\frac{d}{2}\le u\le \frac{2w+d}{2},\hfill \\ 0,\hfill & \mathrm{other}.\hfill \end{array}\right.$$

Here, w denotes the slit width, and d represents the slit distance, both of which are set to $210\mathsf{\mu }$m and $420\mathsf{\mu }$m, respectively. Additionally, ${\rho }_i$ signifies the amplitude of the $i_{th}$ Gaussian plate slit, while ${\tau }_i$ denotes the width parameter governing the phase distribution for the $i_{th}$ Gaussian plate slit.

The parameters associated with the RDSGPP, as depicted in Figure 3, are ${\rho }_1=1,$ ${\rho }_2=0.5, {\tau }_1=100\mathsf{\mu }\mathrm{m},{\tau }_2=65\mathsf{\mu }\mathrm{m}$. In the figure, it is noticeable that when the incident angle ${\theta }_i$ is held constant, particularly in cases such as ${\theta }_i=\pi /4$ and ${\theta }_i=3\pi /8$, the retrieval results exhibit only minor changes as the reflective angle ${\theta }_o$ varies between $\pi /8$, $\pi /4$, and $3\pi /8$, respectively. This observation suggests that the impact of the incident angle ${\theta }_i$ is considerably more significant than that of the reflective angle ${\theta }_o$. Consequently, careful consideration is warranted when selecting the value of the incident angle ${\theta }_i$. Although we have simulated only three incident angle cases, the results affirm that satisfactory retrieval outcomes are achieved when ${\theta }_i=\pi /4$.

Moreover, we note that when ${\theta }_i=\pi /4$, the reconstruction result obtained with ${\theta }_o=\pi /4({\theta }_o=\pi /8)$ surpasses that with ${\theta }_o=3\pi /8$ in Figure 3. Consequently, in subsequent simulations, we fix ${\theta }_i={\theta }_o=\pi /4$. Additionally, we simulate two other instances of RDSGPPs to further substantiate the reliability of the RGI phase retrieval method. The reconstructed phases and amplitudes of these plates are presented in Figure 4. The relative difference between the reconstructed phases and the initial phases is also close to zero, affirming the validity of our RGI phase retrieval scheme.

The second example depicted in Figure 5 is the RDS with the reflection coefficient

$$r_2\left(u\right)=\left\{\begin{array}{ccc}\hfill {\psi }_1{e}^{i{\theta }_{{\psi }_1}}& & -\frac{2w+d}{2}\le u\le -\frac{d}{2}\hfill \\ \hfill {\psi }_1{e}^{i{\theta }_{{\psi }_1}}& & \frac{d}{2}\le u\le \frac{2w+d}{2}\hfill \\ \hfill 0& & \mathrm{other}.\hfill \end{array}\right.$$

Here, ${\psi }_i$ and ${\theta }_{{\psi }_i}$ represent the amplitude and phase of the $i_{th}$ slit, respectively. To demonstrate the reliability of our RGI phase retrieval method, we selected three different RDSs, and their reconstructed phases and amplitudes are presented in Figure 5. In Figure 5, it is evident that the solid blue curves, representing the original phase and amplitude, closely align with the dashed red curves, representing the reconstructed phase and amplitude. This alignment indicates a substantial consistency between the reconstructed and original phases (amplitudes), providing additional evidence for the reliability of our phase reconstruction approach.

Finally, we would like to clarify why quantitative evaluation metrics, such as Signal-to-Noise Ratio (SNR), were not used in assessing the reconstruction results. The evaluation of SNR traditionally leans towards emphasizing amplitude information, making it a commonly used metric. SNR is conventionally computed by dividing the average signal value by the standard deviation of noise, with amplitude playing a pivotal role in this calculation. Phase information, on the other hand, is typically utilized to describe signal variations, periodicity, and relative positions, rather than directly influencing the signal’s intensity. Consequently, direct incorporation of phase information into SNR calculations is less prevalent, as phase is not inherently associated with the concept of ‘noise’.

It’s important to underscore that our amplitude reconstruction relies entirely on the established framework of the most traditional ghost imaging scheme. Additionally, we want to draw attention to our observation that the double integration function in MATLAB yielded suboptimal results when applied to integrate Equation (11a). This observation is manifested in the curve of $I_0$ depicted in Figure 2 (where $I_0$ is derived under the condition of $f\left(x_1\right)=1$, representing the outcome of the most traditional ghost imaging system). According to traditional ghost imaging theory, $I_0$ can be regarded as the square of the object’s amplitude. Ideally, if the double integration function in MATLAB were perfect, the curve of $I_0$ for each slit in Figure 2 should be parallel to the axis. However, it is evident that the curves of $I_0$ for each slit in Figure 2 do not align with the axis; instead, they exhibit significant fluctuations. This discrepancy indicates that the double integration function in MATLAB is not an ideal fit for simulating Equation (11a). Consequently, the reconstructed amplitude based on this foundation is also imperfect (e.g., the reconstructed amplitude in Figure 4 does not perfectly align with the original amplitude information). Similarly, we can infer that $I_1$, $I_2$, $I_3$, $I_4$, and $I_5$ obtained using this double integration function are not flawless, resulting in a deviation between the reconstructed phase information and theoretical values (as evident in simulation results, such as Figure 4).

Despite these challenges, it is crucial to highlight that even though $I_0$, $I_1$, $I_2$, $I_3$, $I_4$, and $I_5$ obtained using the MATLAB double integration function are not perfect, and the simulation results (Figure 4 and Figure 5) maintain a high level of persuasiveness. This resilience indicates that the FSPS method remains a reliable approach.

4. Conclusions

In conclusion, through theoretical derivation and extensive simulations, we have demonstrated that the FSPS method can be directly applied to RGI for obtaining the phase distribution of reflecting objects. This method successfully reconstructed the reflection phases and amplitudes of three complex RDSGPPs and three complicated RDSs. We have explored different combinations of incident angles ${\theta }_i$ and reflective angles ${\theta }_o$ in simulations to analyze their effects on phase reconstruction performance. The simulation results reveal that the incident angle has a significantly greater impact than the reflective angle on phase retrieval results, emphasizing the importance of carefully choosing the incident angle. Our scheme exhibits robust performance when the incident angle is set to ${\theta }_i=\pi /4$. Additionally, under the conditions of ${\theta }_i=\pi /4$, the phase reconstruction result with ${\theta }_o=3\pi /8$ is less favorable than with ${\theta }_o=\pi /4$ (${\theta }_o=\pi /8$), suggesting a preference for smaller reflective angles (${\theta }_o\le \pi /4$).

Crucially, the FSPS method sets itself apart from commonly used GI phase retrieval methods by eliminating the necessity for entangled photon pairs, avoiding the need for complex optical systems, skipping any Fourier transform steps, and foregoing lengthy iterative processes in its phase retrieval procedure. This distinction marks a departure from conventional GI phase retrieval methods. Furthermore, the FSPS method proves successful in acquiring object phase information, whether in TGI [34] or RGI, underscoring its remarkable versatility. We posit that this method is not only applicable to GI but also harbors substantial potential in various other phase retrieval applications.