A phase recovery method

Abstract

This invention belongs to the field of phase recovery technology in optical measurement. It proposes a phase recovery method based on multiple weighted iterations of the light intensity transmission equation using angular spectrum theory. Two light intensity images are collected during light propagation. The light intensity transmission equation is solved using a fast Fourier transform. The result of the light intensity transmission equation is used as the initial phase. A weighted iterative algorithm is then used to perform multiple iterations to invert the measured wavefront phase. This method solves the problems of slow convergence and low accuracy of traditional iterative algorithms. The recovered phase image has smaller errors and significantly improves convergence speed compared to traditional phase recovery algorithms.

Description

Technical Field

[0001] The present invention belongs to the technical field of phase recovery in optical measurement, and specifically designs a phase recovery method of a light intensity transmission equation based on multiple weighted iterations of angular spectrum theory. Background Art

[0002] Light information is composed of phase and intensity. Currently, light intensity detectors such as CCD can directly measure the intensity distribution of light, but cannot measure the phase distribution of light. In order to obtain all the information of the light field, interferometric measurement methods such as holography and shearing interferometry have been proposed. Due to the defects of interferometric measurement methods such as high environmental requirements and the need to add additional reference light, phase retrieval algorithms have been proposed. Since phase recovery is a non-interferometric quantitative measurement, it does not require the introduction of reference light. The measurement optical path structure is relatively simple, and the environmental requirements are lower than those of interferometric measurement, and it has stronger adaptability to the environment. Therefore, phase retrieval algorithms have been widely used in surface profiling, holographic technology, binary optical element design, image encryption, aberration correction of microscopic imaging, large-aperture wavefront detection, etc.

[0003] The main methods of phase retrieval technology are iteration and the transmission intensity equation (TIE). The iteration method, first proposed by R.W. Gerchberg and W.O. Saxton to solve electron microscopy imaging problems, is the Gerchberg-Saxton (GS) algorithm. The iteration method uses measured data to constrain the spatial and frequency domains so that the phase gradually approaches the correct value. The classic GS iteration method simultaneously measures the light intensity distribution at the front and rear focal planes of the lens and repeatedly replaces the amplitude value during the iteration. The GS algorithm, based on angular spectrum theory, iterates based on the diffraction transformation of two adjacent light intensity distributions in the Fresnel diffraction field.

[0004] The TIE method was derived by MRTeague et al. using the Helmholtz equation and the Green's function to provide a solution for TIE. TIE uses the intensity distribution change in the propagation direction of the light field to derive the phase distribution in the perpendicular direction. Currently, there are several main solutions for solving TIE:

[0005] (1) Green's function solution based on Derek Ray boundary conditions, proposed by Teague et al. in 1983;

[0006] (2) Zernike polynomial method, proposed by TE Gureyev and KANugent in 1995;

[0007] (3) Multigrid method, proposed by LJ Allen and MPOxley in 2001;

[0008] (4) Fourier transform method, first applied by Kazuichi Ichikawa to solve TIE.

[0009] The GS iteration method based on angular spectrum theory is an error reduction algorithm, but its convergence speed is too slow and requires a lot of computing resources when used; the numerical accuracy and morphological characteristics of the calculation results of the light intensity transmission equation method are not satisfactory. Summary of the Invention

[0010] To address the shortcomings of iterative and light intensity transfer equation methods in phase retrieval technology, this paper proposes an algorithm (GS-TIE) based on the angular spectrum weighted iteration (GS) algorithm combined with the light intensity transfer equation (TIE). This algorithm quickly and efficiently solves the problems of insufficient accuracy and slow convergence speed when these two algorithms are used alone, and can achieve efficient and accurate phase retrieval. To achieve the above objectives, the present invention adopts the following specific technical solutions:

[0011] A phase recovery method, comprising:

[0012] S1: Use a semi-transparent and semi-reflective mirror to split the light to be measured into two beams, and use the image sensor to collect the light intensity distribution I of the two beams g , I k , and I g with I k The optical path difference between z ;

[0013] S2: The two light intensity distributions I g , I k , the Fourier transform method is used to solve the light intensity transmission equation and obtain the phase distribution P;

[0014] S3: Use the phase distribution P as the initial value of the angular spectrum weighted iterative algorithm and use two light intensity distributions I g , I k Perform weighted iteration to invert the accurate phase distribution.

[0015] Preferably, step S3 specifically includes:

[0016] S301: Calculate f using formula (1) p The complex amplitude is obtained by weighting the amplitude part of

[0017]

[0018] Among them, f p is the complex amplitude obtained by GS iteration;

[0019] f obtained in the previous iteration p Phase;

[0020] α and β are weighting coefficients;

[0021] S302: Light Field The Fresnel diffraction method is used to simulate the propagation of light in free space, that is, the propagation of the light to be measured in free space Δ is calculated by angular spectrum theory. z Complex amplitude distribution after distance:

[0022]

[0023] Among them, g p for In free space propagation Δ z The complex amplitude obtained after

[0024] g p Phase;

[0025] S303: Use the square root of light intensity and the |g calculated in S302 p | combination instead of |g p |, the weighting coefficients α and β remain unchanged, and the light field g is obtained p The weighted complex amplitude:

[0026]

[0027] S304: Yes Perform inverse angular spectrum transformation to obtain the light field Back propagation of Δ in free space z The complex amplitude after:

[0028] Determine whether the number of iterations meets the requirements. If the maximum number of iterations is reached, stop the iteration output. If not satisfied, then p=p+1 and S301 is repeated.

[0029] Preferably, the phase distribution P in step S2 is obtained using the following formula:

[0030]

[0031] in, Z is the propagation distance; λ is the wavelength of the light wave;

[0032] k ⊥ 2 =k x 2 +k y 2 ;k x 、k y are the wave vector components in the x and y directions.

[0033] Preferably, the error SSE between each phase distribution and the true phase is calculated using the following formula:

[0034]

[0035] in, Represents the light intensity of the first surface at the pth iteration.

[0036] Preferably, the optimal parameters of the weighting coefficients α and β determined through multiple calculations are α=1.2 and β=0.2.

[0037] The present invention can achieve the following technical effects:

[0038] 1. The present invention simultaneously adopts two commonly used phase recovery algorithms to solve the problems of slow convergence and local minimum of the GS weighted iterative algorithm based on angular spectrum theory, and low precision and instability of the TE method.

[0039] 2. By combining the two algorithms, the recovered phase image has advantages in both accuracy and acquisition speed.

[0040] 3. The present invention only requires the target light intensity distribution and the light intensity distribution at a specific small distance from it. The required experimental device is easy to build, and the inverted phase result has many advantages over traditional algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 This is a schematic diagram of light intensity sampling of a phase recovery method according to an embodiment of the present invention;

[0042] Figure 2 is a flow chart of an embodiment of the present invention;

[0043] Figure 3 is a target phase to be recovered according to an embodiment of the present invention;

[0044] Figure 4 is the target light propagation Δ of one embodiment of the present invention z After the phase;

[0045] Figure 5 This is an embodiment of the present invention, wherein the light intensity transmission algorithm is used to recover the obtained phase;

[0046] Figure 6 The phase is obtained using a weighted iterative algorithm based on angular spectrum theory according to an embodiment of the present invention;

[0047] Figure 7 This is an embodiment of the present invention, wherein the phase is recovered using a combined algorithm;

[0048] Figure 8 This is the error comparison chart of GS algorithm and combined algorithm. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and do not constitute a limitation of the present invention.

[0050] The purpose of the present invention is to provide a fast and efficient method that solves the problems of insufficient precision and slow convergence when using these two algorithms alone, enabling efficient and accurate phase recovery. The following describes a phase recovery method for the light intensity transmission equation based on multiple weighted iterations of angular spectrum theory, as provided by the present invention, in detail through specific examples.

[0051] First, a semi-transparent and semi-reflective mirror is used to obtain an optical path difference of Δ z The two beams of light are used to collect the light intensity distribution I g , I k ; Then the Fourier transform method is used to solve the light intensity transmission equation to obtain its phase distribution P; Finally, the phase distribution P is used as the initial value of the angular spectrum weighted iterative algorithm, and the two light intensity distributions I g , I k Perform weighted iteration to invert the accurate phase distribution.

[0052] like Figure 2 As shown in the flow chart, the light intensity I is obtained. g , I k , and the phase distribution obtained by the TIE algorithm amplitude After the parameters are equal, the first GS iteration is performed to obtain the With known amplitude Combined into complex amplitude That is, formula (1):

[0053] The angular spectrum is transformed to obtain In free space propagation Δ z The complex amplitude obtained is formula (2):

[0054] Use the square root of the light intensity and the |g calculated in the previous step p | combination instead of |g p |, taking the same α and β, we can get the weighted complex amplitude That is, formula (3);

[0055] Perform inverse angular spectrum transformation on equation (3) to obtain Back propagation of Δ in free space zThe complex amplitude after: That is, formula (4);

[0056] Determine whether the number of iterations p is the maximum value at this time. If not, repeat steps (1), (2), (3), and (4).

[0057] If so, stop the iteration and output the phase

[0058] Figure 1 The schematic diagram of sampling the light intensity is shown in FIG. 1 . In a preferred embodiment of the present invention, a beam with a wavelength of λ = 632.8 nm and a phase distribution according to The target light is divided into two optical paths with an optical path difference Δ by the semi-transparent and semi-reflective mirror BS. z =2mm two beams of light, and use two image sensors to collect their light intensity I g , I k The image sensor used is a CMOS image sensor or a CCD image sensor. In a preferred embodiment of the present invention, the acquisition is achieved by cameras CMOS1 and CMOS2. The two light intensities I collected in the object space are g , I k , the Fourier transform method is used to solve the light intensity transmission equation to obtain the phase distribution P.

[0059] In another embodiment of the present invention, the Fresnel diffraction formula for light propagating in free space satisfies:

[0060]

[0061] Where u(x,y) is the complex amplitude, u0(x,y) is the complex amplitude distribution of the light wave after it propagates a distance z in free space, λ is the wavelength of the light wave, and k is the wave vector of the light.

[0062] It can be verified that u strictly satisfies the parabola equation:

[0063]

[0064] Transforming equation (8) yields the light intensity transmission equation:

[0065]

[0066] Where, I represents light intensity; P represents phase distribution;

[0067] Performing Fourier transform on the optical transmission equation (9) yields:

[0068]

[0069] in, Z is the propagation distance; λ is the wavelength of the light wave;

[0070] k ⊥ 2 =k x 2 +k y 2 ;k x 、k y are the wave vector components in the x and y directions.

[0071] A program is written using formula (5), and the algorithm is named TIE method to calculate the phase distribution calculated by the light intensity transmission equation method.

[0072] In a preferred embodiment of the present invention, the weighting coefficients α=1.2 and β=0.2 are determined after multiple calculations;

[0073] When performing the first iteration, set the number of iterations p = 0 and the initial amplitude to Initial phase at this time That is, formula (5);

[0074] That is, in the first iteration,

[0075] Next, use formula (2) to calculate |g at this time p |;

[0076]

[0077] Use the square root of the light intensity and the |g calculated in the previous step 0 | combination instead of |g 0 |, get g 0 Weighted complex amplitude That is, using formula (3) we can get:

[0078]

[0079] Again, yes Perform inverse angular spectrum transformation to obtain Back propagation of Δ in free space z The complex amplitude f 0+1 :

[0080] Select As the phase of the second iteration,

[0081] That is, at this iteration p=1, the amplitude is

[0082] The second iteration is performed using formula (1):

[0083]

[0084] And so on, p=p+1 until the number of iterations meets the requirement, stop iteration, and output

[0085] In another embodiment of the present invention, the default number of iterations is set to 100, the sampling plane is a square with a length of 5 mm, the number of sampling points is 128, the amplitude of the target light is randomly set, and the phase of the target light field is generated by computer simulation, as shown in FIG. Figure 3 As shown;

[0086] According to the angular spectrum propagation theory, the phase distribution of the target light after propagation is obtained, such as Figure 4 As shown;

[0087] According to the solution of the light intensity transmission equation obtained in S2, the phase distribution of the target light is solved using the TIE algorithm. The result is as follows Figure 5 As shown;

[0088] Directly use the GS iteration method based on angular spectrum theory, randomly set the initial phase distribution (that is, do not use the phase obtained by the TIE algorithm as the initial phase), iterate 100 times, and the resulting phase distribution is as follows Figure 6 As shown;

[0089] Using the S3 algorithm, namely the GS-TIE algorithm, the results obtained after 100 iterations are as follows Figure 7 As shown;

[0090] In a preferred embodiment of the present invention, the average value SSE of the difference between the phase obtained by each algorithm and the actual phase is calculated in the program, and the error curve is drawn using formula (6) as shown in FIG. Figure 8 As shown;

[0091] For this phase distribution, it can be clearly seen that the GS-TIE algorithm has obvious advantages over the two algorithms:

[0092] Since the initial phase distribution of the GS-TIE algorithm is provided by the TIE algorithm, Figure 8 The initial point of the GS-TIE algorithm error is the error of the TIE algorithm, which is smaller than that of the GS algorithm. The error of the GS algorithm is much larger than the errors of the TIE and GS-TIE algorithms throughout the entire iteration process. Figure 8 It can be seen that the error decreases significantly in the first 10 iterations of the GS-TIE algorithm, demonstrating that the GS-TIE algorithm is more accurate than the TIE algorithm and recovers a more precise phase. Furthermore, compared to the GS algorithm, the GS-TIE algorithm reduces the error more quickly due to the relatively accurate initial phase distribution, requiring fewer iterations.

[0093] In summary, the GS-TIE algorithm takes into account the advantages of both GS and TIE algorithms, and provides a fast and accurate algorithm for solving the phase recovery problem.

[0094] In the description of this specification, the reference terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" mean that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art can combine and combine different embodiments or examples described in this specification and features of different embodiments or examples without contradiction.

[0095] Although the embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention.

[0096] The above specific embodiments of the present invention do not constitute a limitation on the scope of protection of the present invention. Any other corresponding changes and modifications made based on the technical concept of the present invention should be included in the scope of protection of the claims of the present invention.

Claims

1. A phase recovery method, characterized in that: include: S1: Use a semi-transparent and semi-reflective mirror to split the light to be measured into two beams, and use an image sensor to collect the light intensity distribution of the two beams. 、 ,and and The optical path difference between ; S2: The two light intensity distributions 、 , the Fourier transform method is used to solve the light intensity transmission equation and obtain the phase distribution P The phase distribution P in step S2 is obtained using the following formula: （5） in, ; Z is the propagation distance; λ represents the wavelength of the light wave; ; k x 、 k y are the wave vector components in the x and y directions; S3: The phase distribution P As the initial value of the angular spectrum weighted iterative algorithm, and using the two light intensity distributions 、 Perform weighted iteration to invert the accurate phase distribution.

2. The phase recovery method according to claim 1, characterized in that: Step S3 specifically includes: S301: Use formula (1) to calculate The complex amplitude is obtained by weighting the amplitude part of ; （1） in, is the complex amplitude obtained by GS iteration; The result of the previous iteration Phase; α 、 β is the weighting coefficient; S302: For the complex amplitude The Fresnel diffraction method is used to simulate the propagation of light in free space, that is, the propagation of the light to be measured in free space is calculated by the angular spectrum theory. Complex amplitude distribution after distance: （2） in, g p For the Propagation in free space The complex amplitude obtained after for g p Phase; S303: Use the square root of light intensity and the value calculated in S302 Combination of , the weighting coefficient α 、 β unchanged, and the complex amplitude is obtained g p The weighted complex amplitude: （3） S304: Yes Perform inverse angular spectrum transformation to obtain the complex amplitude Back propagation in free space The complex amplitude after: （4） Determine whether the number of iterations meets the requirements. If the maximum number of iterations is reached, stop the iteration output. If not satisfied, , repeat S301.

3. The phase recovery method according to claim 1, wherein: The error SSE between each phase distribution and the true phase is calculated using the following formula: （6） in, Representative p The light intensity of the first face at iteration .

4. The phase recovery method according to claim 2, wherein: The weighting coefficient α 、 β The optimal parameters determined after multiple calculations are α =1.2, β =0.2.

Images