A multi-azimuth differential phase reconstruction method

Abstract

This invention discloses a multi-directional differential phase reconstruction method. For synchronous phase-shifted transverse shearing interferograms in two directions, a phase-correlated image registration algorithm is used to match the positions of the shearing interferograms in both directions. A four-step phase-shifting algorithm is used to obtain the wrapped phase in both directions. An optimized DCT-based least-squares phase unwrapping algorithm is used to unwrap the wrapped phase. The shearing amount of the shearing interferogram is calculated through edge detection and fitting. The wavefront to be measured is reconstructed using a least-squares method based on differential Zernike polynomials. An anti-tilt algorithm eliminates the introduced tilt error, completing the differential phase reconstruction. This invention utilizes an optimized unwrapping algorithm to overcome the unwrapping failure problem caused by noise, making the unfolded phase smoother, more continuous, and closer to the true phase. Furthermore, data is processed in the circular domain, eliminating the need for Zernike polynomial orthogonalization. This invention can be applied to processing multi-directional transverse shearing interferometry.

Description

Technical Field

[0001] This invention belongs to the field of wavefront detection technology, specifically a multi-directional differential phase reconstruction method for detecting transverse shear interference testing technology of the surface shape of transmitted wavefront and optical element. Background Technology

[0002] In optical systems, the quality of the surface shape of optical elements directly affects the imaging quality of the system. High-precision wavefront detection technology can improve the manufacturing quality of optical elements, thereby ensuring the imaging quality of the entire optical system. In the field of optical testing, wavefront detection plays an important role. By testing the wavefront of optical elements, the surface shape information, optical uniformity, and wavefront aberrations of the optical elements or system can be directly obtained from the obtained wavefront information. Based on these parameters, the optical quality of the optical elements or system can be evaluated.

[0003] In the precision polishing stage of optical component manufacturing, optical interferometry is often used to detect surface shape and transmitted wavefront. Shearing interferometry is a measurement technique that utilizes the interference between the intrinsic light wave and its replica light wave. Since the intrinsic and replica light waves are spatially misaligned, no reference light wave is needed, resulting in a simple optical path structure. This eliminates the influence of the reference mirror's surface shape error on the detection results. Furthermore, shearing interferometry systems generally use common-path phase shift measurements, eliminating the influence of environmental vibrations and air disturbances on the interferometric measurement. Combined with synchronous phase shift technology, absolute common-path phase shift measurement of the measured phase can be achieved, eliminating the influence of environmental vibrations and air disturbances on the interferometric measurement. However, due to the self-comparison characteristic of transverse shearing interferometry wavefronts, the wavefront being measured forms an interferogram of the differential wavefront after passing through the shearing interferometer. After phase unwrapping, the resulting image is the phase distribution of the differential wavefront in the shearing direction of the measured wavefront, not the measured wavefront distribution carrying the surface shape of the measured component. Therefore, wavefront reconstruction using a specific algorithm is required to obtain the phase distribution of the measured surface. Therefore, research on multi-directional differential phase reconstruction technology for synchronous phase-shifting shearing interference is of great significance for the detection of the surface shape and transmitted wavefront of optical elements.

[0004] Current wavefront reconstruction algorithms mainly include the region method and the mode method. The region method obtains the one-dimensional distribution of the wavefront to be measured by summing and inversely operating on the differential phase values, and then uses the least squares method to fit the two-dimensional distribution of the wavefront phase. It has the advantages of simple algorithm and fast operation, but it requires multiple measurements, which complicates the structure of the shearing interferometer and the measurement process. At the same time, the measurement repeatability and measurement accuracy will also decrease. The mode method expands the wavefront phase within the full aperture into different modes, and then uses the measurement data within the full aperture to fit the optimal expansion coefficients, thereby obtaining the complete wavefront expansion, and finally reconstructing the wavefront to be measured. It has the advantages of strong noise resistance and high computational efficiency, but it omits higher-order terms and loses high-frequency information of the wavefront being measured. The expansions include Legendre polynomials, power polynomials, Fourier series, Zernike polynomials, etc. Among them, Zernike polynomials are often used in wavefront fitting due to their good convergence, high fitting accuracy of optical wavefronts, short calculation time, shearing amount not limited by sampling interval, and easy connection with aberrations. Summary of the Invention

[0005] The purpose of this invention is to address the problems existing in the prior art by providing a multi-directional differential phase reconstruction method to meet the needs of synchronous phase-shifting shearing interferometry testing technology. This method can obtain the phase distribution to be measured from multi-directional phase-shifting shearing interferograms, achieve high-precision measurement of the surface shape of the transmitted wavefront and optical elements, and can be well applied to wavefront dynamic testing.

[0006] The technical solution to achieve the objective of this invention is as follows: On one hand, a multi-directional differential phase reconstruction method is provided, the method comprising the following steps:

[0007] Step 1: Register the phase-shifted transverse shearing interferograms in the x and y directions respectively;

[0008] Step 2: Obtain the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm;

[0009] Step 3: Transform the phase expansion into the optimal solution of the function, and use the least squares method of weighted iterative DCT to unwrap the differential phase and obtain the differential phase distribution in the x and y directions;

[0010] Step 4: Perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the center of the circle and take the average value to obtain the shearing amount;

[0011] Step 5: Based on the differential phase distributions in the x and y directions obtained in Step 3 and the shearing amount obtained in Step 4, wavefront reconstruction is performed using the least squares method based on the differential Zernike polynomial.

[0012] Further, the position registration of the phase-shifted transverse shearing interferograms in the x and y directions, as described in step 1, specifically includes:

[0013] Step 1-1: Obtain the translation amount of each phase-shifted transverse shearing interferogram relative to the first phase-shifted transverse shearing interferogram using a phase-correlation image registration algorithm;

[0014] Steps 1-2: Based on the translation amount, perform reverse translation on the phase-shifting transverse shearing interferograms to be matched so that the positions of each phase-shifting transverse shearing interferogram are matched.

[0015] Further, the phase-correlation-based image registration algorithm described in step 1-1 obtains the translation amount of each phase-shifting transverse shearing interferogram relative to the first phase-shifting transverse shearing interferogram, specifically including:

[0016] Let f1 be the reference image of the phase-shifted transverse shearing interferogram, and f2 be the image to be registered that has been translated in the spatial domain. Then:

[0017] f1 = f(x, y)

[0018] f2 = f(x - Δx, y - Δy)

[0019] In the formula, (Δx, Δy) represents the translation of the image to be registered f2 relative to the reference image f1 along the x-axis and y-axis, and f(x, y) represents the pixel value of the image at (x, y).

[0020] Let the size of each phase-shifted transverse shearing interferogram be M×N. Performing Fourier transforms on f1 and f2 respectively, we have:

[0021] F1(u, v) = F1(u, v)

[0022]

[0023] In the formula, F1(u,v) and F1(u,v) are the spectrum diagrams after the Fourier transform of f1 and f2, respectively;

[0024] Calculate the normalized power spectra of the two in the frequency domain and perform an inverse Fourier transform to obtain the unit impulse function C(x, y):

[0025]

[0026] In the formula, * represents the conjugate operation, and δ is the Dirac function;

[0027] The translation amount can be obtained by determining the coordinates corresponding to the maximum value of the unit impulse function:

[0028]

[0029]

[0030] Furthermore, step 2, which involves obtaining the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm, specifically includes:

[0031] I1=A+Bcos(φ+π / 2)=A-Bsin(φ)

[0032] I2=A+Bcos(φ+π)=A-Bcos(φ)

[0033] I3=A+Bcos(φ+3π / 2)=A+Bsin(φ)

[0034] I4=A+Bcos(φ+2π)=A+Bcos(φ)

[0035]

[0036] In the formula, I1, I2, I3, and I4 are the light intensities of the interference field, A is the background light intensity of the interference field, and B is the modulation index of the interference field.

[0037] Further, step 3 involves transforming the phase expansion into solving for the optimal solution of the function, and using the least squares method of weighted iterative DCT to unwrap the differential phase, obtaining the differential phase distributions in the x and y directions. The specific process includes:

[0038] Step 3-1: Represent the least squares problem as a vector form determined by multiple factors:

[0039] Ax = b

[0040] Further expressed using cosine transform:

[0041] A T Ax = A T b

[0042] Where x is the phase value length MN result vector, b is the wrapped phase difference vector of length N(M-1)+M(N-1), and T represents matrix transpose;

[0043] Step 3-2, the weighted least squares problem is expressed as:

[0044] WAx = Wb

[0045] Further expressed using the discrete cosine transform, it can be represented as:

[0046] A T W T WAx = A T W T Wb

[0047] Where W is the weight matrix associated with the pixel weights;

[0048] Let Q = A T W T WA,

[0049] Then we have:

[0050] Qφ=c

[0051] Where b is the weighted phase difference vector, and c is the weighted discrete phase Laplace vector of the phase difference.

[0052] Decompose matrix Q into matrix P and difference matrix D;

[0053] Step 3-3: Calculate c using the initial weight data and the wrapped phase difference vector;

[0054] Steps 3-4: Set the initial conditions as iteration number k = 0, differential phase φ k =0, maximum number of iterations k max ;

[0055] Steps 3-5: Iteratively calculate the differential phase φ k+1 :

[0056] φ k+1 =c-Dφ k

[0057] Steps 3-6: Solve for ρ using the DCT least squares method. k :

[0058] Pφ k+1 =ρ k

[0059] In the formula, P is a vector φ k+1 The matrix to which the discrete Laplace operation is performed, ρ k It is a vector that contains a discrete Laplace operation on the phase difference of the package;

[0060] Steps 3-7: Determine if k < k max If the condition is true, continue iterating from step 3-3 to step 3-4 until k = k. max Finally, we can solve for φ. k+1 .

[0061] Further, step 4 involves edge detection and fitting of the phase-shifted transverse shearing interferograms in the x and y directions, calculating the distance between the center points, and taking the average value to obtain the shearing amount. Specifically, this includes:

[0062] Step 4-1: Perform edge detection on the shearing interferogram, and then fill the gaps between the stripes;

[0063] Step 4-2, perform edge detection again;

[0064] Step 4-3: Fit two spot circles of the sheared beam from the left and right edges of the phase-shifted transverse shearing interference pattern in the x direction or the upper and lower edges in the y direction, respectively. Calculate the distance between the centers of the two circles and take the average value to obtain the magnitude of the shearing.

[0065] Furthermore, in step 5, based on the differential phase distributions in the x and y directions obtained in step 3 and the shearing amount obtained in step 4, wavefront reconstruction is performed using the least squares method based on the difference Zernike polynomial. The specific process includes:

[0066] According to the Zernike polynomial, the wavefront is represented as W(x, y):

[0067] W(x, y) = ∑a i Z i (x, y)

[0068] Among them, Z i (x, y) represents the i-th Zernike polynomial, a i This indicates its corresponding coefficient;

[0069] The shear wavefronts along the x and y directions are represented by ΔW, respectively. x (x, y) and ΔW y (x, y):

[0070] ΔW x (x,y)=W(x+s,y)-W(x,y)=∑a xi [Z i (x+s,y)-Z i [x, y]=∑a xi Z i (x, y)

[0071] ΔW y (x,y)=W(x,y+s)-W(x,y)=∑a yi [Z i (x, y+s)-Z i [x, y]=∑a yi Z i (x, y)

[0072] Where s is the shear rate;

[0073] The above relationship can be expressed in vector form:

[0074]

[0075] Right now

[0076] In the formula, a x =[a x1 a x2 , ..., a xi ] T a y =[a y1 a y2 , ..., a yi ] T a x1 and a y1 Let T represent the Zernike coefficients of the i-th shear wavefront, respectively. x and T y The coefficient transformation matrix is ​​obtained by combining the relationship between the difference Zernike polynomials and the Zernike polynomials and eliminating the correlation columns.

[0077]

[0078] The Zernike coefficients before the measured wavefront are expressed as:

[0079]

[0080] Among them, T + It is the generalized inverse matrix of T, T + =(T + T - )T T ;

[0081] The wavefront under test is fitted with the Zernike coefficients of the wavefront under test, thus completing the reconstruction from the shear wavefront to the wavefront under test.

[0082] On the other hand, a multi-directional differential phase reconstruction system is provided, the system comprising:

[0083] The first module is used to perform position registration of the phase-shifted transverse shearing interferograms in the x and y directions, respectively.

[0084] The second module is used to obtain the differential phases wrapped in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm.

[0085] The third module is used to transform the phase expansion into the optimal solution of the solution function. Based on the least squares method of weighted iterative DCT, the differential phase is unwrapped to obtain the differential phase distribution in the x and y directions.

[0086] The fourth module is used to perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the centers of the circles and take the average value to obtain the shearing amount;

[0087] The fifth module is used to perform wavefront reconstruction based on the differential phase distribution and shearing amount in the x and y directions obtained above, using the least squares method based on the differential Zernike polynomial.

[0088] Compared with the prior art, the significant advantages of this invention are:

[0089] (1) An optimized least squares unpacking algorithm based on weighted iterative DCT is proposed, which can process noise and residual information more smoothly. When processing low-frequency information, sin(φ) = φ is used to further suppress the unpacking failure caused by noise.

[0090] (2) The Canny operator is used to detect and segment the edges of the shearing interferogram to calculate the shearing amount. The average value is obtained by processing the shearing interferogram from multiple directions, which can be accurate to the pixel level and improve the accuracy of obtaining the shearing amount.

[0091] (3) The least squares method based on the differential Zernike polynomial is used for wavefront reconstruction. The shearing interferogram is processed within the circular mask area. Therefore, the Zernike polynomial is still orthogonal and does not need to be Gram-Schmidt orthogonalized. It is not limited by the shearing amount or sampling interval. The missing high-frequency information can be almost ignored by fitting the 36 Zernike polynomials, which can achieve high-precision differential phase reconstruction.

[0092] The present invention will now be described in further detail with reference to the accompanying drawings. Attached Figure Description

[0093] Figure 1 This is a flowchart of the multi-directional differential phase reconstruction method of the present invention.

[0094] Figure 2 Here is an example of a simulation of a synchronous phase-shifting shearing interferogram in two directions, wherein... Figure 2 (a) in the diagram represents four phase-shifting shearing interferograms in the x-direction. Figure 2 (b) in the diagram represents four phase-shifted shearing interferograms in the y-direction.

[0095] Figure 3 This is a sheared interferogram after position registration using a phase-correlation-based image registration algorithm in one embodiment. Figure 3 (a) shows four shearing interferograms after registration in the x-direction, where Figure 3 (b) in the figure shows four shearing interferograms after registration in the y direction.

[0096] Figure 4 In one embodiment, the differential phase encompassing two directions is obtained using a four-step phase-shifting algorithm, wherein... Figure 4In the figure, (a) represents the differential phase enclosed in the x-direction. Figure 4 In the equation (b), the differential phase enclosed in the y-direction is represented.

[0097] Figure 5 In one embodiment, the differential phases of the two directions are obtained by unwrapping using the least squares method based on weighted iterative DCT. Figure 5 In the figure, (a) represents the differential phase after unwrapping in the x-direction. Figure 5 In the equation (b), the differential phase after unwrapping in the y direction is represented.

[0098] Figure 6 This is a schematic diagram illustrating the process of obtaining shearing amount using an edge detection and fitting algorithm in one embodiment.

[0099] Figure 7 This section compares the final reconstructed phase simulation results and reconstruction errors using the method of the present invention in one embodiment. Figure 7 (a) in the image represents the final reconstructed phase result. Figure 7 (b) in the figure represents the simulation reconstruction error.

[0100] Figure 8 This is the final reconstructed phase of a real phase-shifting shearing interferogram processed using the method of the present invention in one embodiment. Figure 8 (a) in the figure is the phase-shifting shearing interference pattern obtained in the x-direction from the experiment. Figure 8 (b) in the figure is the phase-shifting shearing interferogram obtained in the y-direction from the experiment. Figure 8 (c) in the image is the final reconstructed 3D wavefront plot. Figure 8 (d) in the diagram represents the final reconstructed two-dimensional wavefront plot. Detailed Implementation

[0101] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0102] It should be noted that if the embodiments of the present invention involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of the components in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0103] Furthermore, if the embodiments of this invention involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the technical solutions of the various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by this invention.

[0104] In one embodiment, combined Figure 1 A multi-directional differential phase reconstruction method is provided, the method comprising the following steps:

[0105] Step 1: Register the phase-shifted transverse shearing interferograms in the x and y directions respectively;

[0106] Step 2: Obtain the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm;

[0107] Step 3: Transform the phase expansion into the optimal solution of the function, and use the least squares method of weighted iterative DCT to unwrap the differential phase and obtain the differential phase distribution in the x and y directions;

[0108] Step 4: Perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the center of the circle and take the average value to obtain the shearing amount;

[0109] Step 5: Based on the differential phase distributions in the x and y directions obtained in Step 3 and the shearing amount obtained in Step 4, wavefront reconstruction is performed using the least squares method based on the differential Zernike polynomial.

[0110] Furthermore, in one embodiment, the phase-shifted transverse shearing interferograms in the x and y directions described in step 1 are obtained by using a phase-correlated image registration algorithm to obtain the translation amount of each interferogram relative to the first interferogram, and then performing an inverse translation to match the positions of each shearing interferogram. The specific process is as follows:

[0111] Step 1-1: Obtain the translation amount of each phase-shifted transverse shearing interferogram relative to the first phase-shifted transverse shearing interferogram using a phase-correlation image registration algorithm; specifically including:

[0112] Let f1 be the reference image of the phase-shifted transverse shearing interferogram, and f2 be the image to be registered that has been translated in the spatial domain. Then:

[0113] f1 = f(x, y)

[0114] f2 = f(x - Δx, y - Δy)

[0115] In the formula, (Δx, Δy) represents the translation of the image to be registered f2 relative to the reference image f1 along the x-axis and y-axis, and f(x, y) represents the pixel value of the image at (x, y).

[0116] Let the size of each phase-shifted transverse shearing interferogram be M×N. Performing Fourier transforms on f1 and f2 respectively, we have:

[0117] F1(u, v) = F1(u, v)

[0118]

[0119] In the formula, F1(u,v) and F1(u,v) are the spectrum diagrams after the Fourier transform of f1 and f2, respectively;

[0120] Calculate the normalized power spectra of the two in the frequency domain and perform an inverse Fourier transform to obtain the unit impulse function C(x, y):

[0121]

[0122] In the formula, * represents the conjugate operation, and δ is the Dirac function;

[0123] The translation amount can be obtained by determining the coordinates corresponding to the maximum value of the unit impulse function:

[0124]

[0125]

[0126] Steps 1-2: Based on the translation amount, perform reverse translation on the phase-shifting transverse shearing interferograms to be matched so that the positions of each phase-shifting transverse shearing interferogram are matched.

[0127] Furthermore, in one embodiment, step 2, which involves obtaining the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm, specifically includes:

[0128] I1=A+Bcos(φ+π / 2)=A-Bsin(φ)

[0129] I2=A+Bcos(φ+π)=A-Bcos(φ)

[0130] I3=A+Bcos(φ+3π / 2)=A+Bsin(φ)

[0131] I4=A+Bcos(φ+2π)=A+Bcos(φ)

[0132]

[0133] In the formula, I1, I2, I3, and I4 are the light intensities of the interference field, A is the background light intensity of the interference field, and B is the modulation index of the interference field.

[0134] Furthermore, in one embodiment, step 3 involves transforming the phase expansion into the optimal solution of the solution function, and using the least squares method of weighted iterative DCT to unwrap the differential phase, thereby obtaining the differential phase distributions in the x and y directions.

[0135] Here, the specific approach is as follows:

[0136] Let Ψ be the phase of the discrete points in an M×N rectangular region. i，j , Φ i，j For the corresponding unpacking phase, we have:

[0137] Φ i，j =Ψ i，j +2πk

[0138] Where k is the interference order, and -π≤Ψ i，j ≤π, i=0, 1,..., M-1; j=0, 1,...N-1

[0139] Define the wrapper operator W r We can obtain:

[0140] W r {Φ i，j}=Ψ i，j

[0141] definition:

[0142]

[0143]

[0144]

[0145] Solving the system of equations S in the least squares sense yields the unpacking phase Φ. i，j .

[0146]

[0147] By performing a simple transformation on the normal equations of the least squares system, we obtain the discrete Poisson equations on an M×N rectangular grid:

[0148] Φ i+1，j -2Φ i，j +Φ i-1，j +Φ i，j+1-2Φ i，j +Φ i，j-1 =ρ i，j

[0149] in,

[0150]

[0151] This formula is valid for all rectangular discrete points i = 0, 1, ..., M-1; j = 0, 1, ..., N-1, and has been used to calculate ρ. i，j The phase difference is non-zero only within the rectangular region. The Poisson equation has Neumann boundary conditions.

[0152] If Φ has already been obtained i，j DCT spectrum From the inverse discrete cosine transform (IDCT), we can obtain:

[0153]

[0154] Where ω1(m)=1 / 2, m=0, ω1(m)=1, 1≤m≤M-1

[0155] ω2(n)=1 / 2,n=0 ω1(n)=1,1≤n≤N-1

[0156] After solving and simplifying the equations, we obtain the exact solution in the DCT domain:

[0157]

[0158] Based on this, the algorithm is optimized using a weighted sum and iterative approach to obtain the optimal solution. The aforementioned unweighted least squares algorithm can be represented in vector form under multi-factor determination:

[0159] Ax = b

[0160] Using the cosine transform in the least squares sense, the least squares solution is the solution to the standard equation, which can be expressed as:

[0161] A T Ax = A T b

[0162] Where x is the phase value length MN result vector, b is the wrapped phase difference vector of length N(M-1)+M(N-1), and T denotes matrix transpose. Therefore, the above equation can be rewritten as:

[0163] Pφ=ρ

[0164] Where P = A T A is a vector The matrix for which the discrete Laplace operation is performed, ρ = A Tb is a vector containing the discrete Laplace operation on the phase difference of the encapsulated data. Therefore, the above equation can be solved by performing phase unpacking using the DCT algorithm.

[0165] The weighted least squares problem is represented as:

[0166] WAx = Wb

[0167] Where W is the weight matrix, finding its least squares solution is equivalent to finding the orthogonal equation:

[0168] A T W T WAx = A T W T Wb

[0169] Let Q = A T W T WA, Then there is

[0170] Qφ=c

[0171] in, is the weighted phase difference vector, and c is the discrete phase Laplace vector with weighted phase differences. This equation defines the matrix equation for the 2D weighted least squares phase unpacking problem. The following section provides an algorithm to solve this equation.

[0172] The direct method can be directly applied to solve the unweighted least squares problem, but it cannot be directly applied to solve the weighted least squares problem. An iterative algorithm must be constructed. This scheme is based on a simple Picard iterative algorithm to solve the weighted least squares problem.

[0173] Decompose matrix Q into matrix P and difference matrix D:

[0174] Q = P + D

[0175] Further results were obtained:

[0176] (P+D)φ=c

[0177] or

[0178] Pφ=c-Dφ

[0179] By solving the equation through multiple iterations, we have:

[0180] Pφ k+1 =c-Dφ k =ρ k

[0181] Where k is the number of iterations, and the function is called as follows:

[0182] Dφ k =(QP)φk

[0183] Where Q = A T W T WA is a discrete Laplace operator that corrects for the weighted phase difference, P = A T A is the discrete Laplace operator operating on the unweighted phase difference. The above equation is the most fundamental iterative method for solving the weighted least squares unwrapping problem. ρ k This can be represented as a two-dimensional array:

[0184] ρ i，j =c i，j -[U i，j (φ i+1，j -φ i，j )-U i-1，j (φ i，j -φ i-1，j )+V i，j (φ i，j+1 -φ i，j )-V i，j-1 (φ i，j -φ i，j-1 )]

[0185] Among them, c i，j For the weighted phase Laplace operator.

[0186] In summary, the specific process of step 3 includes:

[0187] Step 3-1: Represent the least squares problem as a vector form determined by multiple factors:

[0188] Ax = b

[0189] Further expressed using cosine transform:

[0190] A T Ax = A T b

[0191] Where x is the phase value length MN result vector, b is the wrapped phase difference vector of length N(M-1)+M(N-1), and T represents matrix transpose;

[0192] Step 3-2, the weighted least squares problem is expressed as:

[0193] WAx = Wb

[0194] Further expressed using the discrete cosine transform, it can be represented as:

[0195] A T W T WAx = A T W T Wb

[0196] Where W is the weight matrix associated with the pixel weights;

[0197] Let Q = A T W T WA,

[0198] Then we have:

[0199] Qφ=c

[0200] in, is the weighted phase difference vector, and c is the weighted discrete phase Laplace vector of the phase difference;

[0201] Decompose matrix Q into matrix P and difference matrix D;

[0202] Step 3-3: Calculate c using the initial weight data and the wrapped phase difference vector;

[0203] Steps 3-4: Set the initial conditions as iteration number k = 0, differential phase φ k =0, maximum number of iterations k max ;

[0204] Steps 3-5: Iteratively calculate the differential phase φ k+1 :

[0205] φ k+1 =c-Dφ k

[0206] Steps 3-6: Solve for ρ using the DCT least squares method. k :

[0207] Pφ k+1 =ρ k

[0208] In the formula, P is a vector φ k+1 The matrix to which the discrete Laplace operation is performed, ρ k It is a vector that contains a discrete Laplace operation on the phase difference of the package;

[0209] Steps 3-7: Determine if k < k max If the condition is true, continue iterating from step 3-5 to step 3-6 until k = k. max Finally, we can solve for φ. k+1 .

[0210] Furthermore, in one embodiment, step 4, which involves edge detection and fitting of the phase-shifted transverse shearing interferograms in the x and y directions, calculating the distance between the centers of the circles, and taking the average value to obtain the shearing amount, specifically includes:

[0211] Step 4-1: Use the Canny operator to perform edge detection on the shearing interferogram, and then fill the gaps between the stripes;

[0212] Step 4-2: Perform edge detection again using the Canny operator;

[0213] Step 4-3: Fit two spot circles of the sheared beam from the left and right edges of the phase-shifted transverse shearing interference pattern in the x direction or the upper and lower edges in the y direction, respectively. Calculate the distance between the centers of the two circles and take the average value to obtain the magnitude of the shearing.

[0214] Furthermore, in one embodiment, step 5, based on the differential phase distributions in the x and y directions obtained in step 3 and the shearing amount obtained in step 4, employs a least-squares method based on differential Zernike polynomials for wavefront reconstruction. The specific process includes:

[0215] According to the Zernike polynomial, the wavefront is represented as W(x, y):

[0216] W(x, y) = ∑a i Z i (x, y)

[0217] Among them, Z i (x, y) represents the i-th Zernike polynomial, a i This indicates its corresponding coefficient;

[0218] The shear wavefronts along the x and y directions are represented by ΔW, respectively. x (x, y) and ΔW y (x, y):

[0219] ΔW x (x,y)=W(x+s,y)-W(x,y)=∑a xi [Z i (x+s,y)-Z i [x, y]=∑a xi Z i (x, y)

[0220] ΔW y (x,y)=W(x,y+s)-W(x,y)=∑a yi [Z i (x, y+s)-Z i [x, y]=∑a yi Z i (x, y) where s is the shear rate;

[0221] The above relationship can be expressed in vector form:

[0222]

[0223] Right now

[0224] In the formula, a x =[a x1 a x2 , ..., a xi ] T a y =[a y1 a y2 , ..., a yi ] T a x1 and a y1 Let T represent the Zernike coefficients of the i-th shear wavefront, respectively. x and T y The coefficient transformation matrix is ​​obtained by combining the relationship between the difference Zernike polynomials and the Zernike polynomials and eliminating the correlation columns.

[0225]

[0226] The Zernike coefficients before the measured wavefront are expressed as:

[0227]

[0228] Among them, T + It is the generalized inverse matrix of T, T + =(T + T - )T T ;

[0229] The wavefront under test is fitted with the Zernike coefficients of the wavefront under test, thus completing the reconstruction from the shear wavefront to the wavefront under test.

[0230] In one embodiment, a multi-directional differential phase reconstruction system is provided, the system comprising:

[0231] The first module is used to perform position registration of the phase-shifted transverse shearing interferograms in the x and y directions, respectively.

[0232] The second module is used to obtain the differential phases wrapped in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm.

[0233] The third module is used to transform the phase expansion into the optimal solution of the solution function. Based on the least squares method of weighted iterative DCT, the differential phase is unwrapped to obtain the differential phase distribution in the x and y directions.

[0234] The fourth module is used to perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the centers of the circles and take the average value to obtain the shearing amount;

[0235] The fifth module is used to perform wavefront reconstruction based on the differential phase distribution and shearing amount in the x and y directions obtained above, using the least squares method based on the differential Zernike polynomial.

[0236] Specific limitations regarding the multi-directional differential phase reconstruction system can be found in the limitations of the multi-directional differential phase reconstruction method described above, and will not be repeated here. Each module in the aforementioned multi-directional differential phase reconstruction system can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in hardware or independently of the processor in a computer device, or stored in software in the memory of a computer device, so that the processor can call and execute the corresponding operations of each module.

[0237] In one embodiment, a computer device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to perform the following steps:

[0238] Step 1: Register the phase-shifted transverse shearing interferograms in the x and y directions respectively;

[0239] Step 2: Obtain the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm;

[0240] Step 3: Transform the phase expansion into the optimal solution of the function, and use the least squares method of weighted iterative DCT to unwrap the differential phase and obtain the differential phase distribution in the x and y directions;

[0241] Step 4: Perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the center of the circle and take the average value to obtain the shearing amount;

[0242] Step 5: Based on the differential phase distributions in the x and y directions obtained in Step 3 and the shearing amount obtained in Step 4, wavefront reconstruction is performed using the least squares method based on the differential Zernike polynomial.

[0243] For specific limitations on each step, please refer to the limitations on the multi-directional differential phase reconstruction method mentioned above, which will not be repeated here.

[0244] As a specific example, in one embodiment, the multi-directional differential phase reconstruction method of the present invention is verified and illustrated.

[0245] In this embodiment, four phase-shifted shearing interferograms in two directions are simulated based on the set wavefront phase distribution, as follows: Figure 2 As shown, phase reconstruction is performed, and the steps are as follows:

[0246] Step 1: For the phase-shifted transverse shearing interferograms in the x and y directions, obtain the translation amount of each interferogram relative to the first interferogram using a phase-correlated image registration algorithm, and then perform the reverse translation as follows: Figure 3 As shown, the positions of each shearing interferogram are matched.

[0247] Step 2: Obtain the differential phase enclosed in the x and y directions from the light intensity of the sheared interferogram using a four-step phase-shifting algorithm, such as... Figure 4 As shown.

[0248] Step 3: Transform the phase expansion into solving for the optimal solution of the function. Based on the least squares method of weighted iterative DCT, unwrap the differential phase to obtain the differential phase distributions in the x and y directions, as shown below. Figure 5 As shown.

[0249] Step 4: Perform edge detection and fitting on a total of 8 interferograms (phase-shifted shearing interferograms in the x and y directions), as follows: Figure 6 As shown, the distance between the centers of the circles is calculated and the average value is taken to obtain the shearing amount.

[0250] Step 5: Based on the differential phases in the x and y directions obtained in Step 3 and the shearing amount obtained in Step 4, the phase distribution obtained by wavefront reconstruction using the least squares method based on the differential Zernike polynomial, and the reconstruction error compared with the simulated phase, are shown below. Figure 7 As shown, the reconstruction error PV and RMS values ​​are very small, with PV being 3.0424e. -3 λ, RMS value is 3.9264e -4 λ.

[0251] The final reconstructed phase obtained by processing the actual phase-shifting shearing interferogram using this method is as follows: Figure 8 As shown.

[0252] In summary, this invention processes the shear interferograms in the x and y directions using a phase-correlation-based image registration algorithm, a four-step phase-shifting algorithm, and a phase unwrapping algorithm based on weighted iterative DCT to obtain the shear phase distribution. Then, it reconstructs the wavefront to be measured using a least-squares wavefront reconstruction algorithm based on differential Zernike polynomials. The results of the embodiments show that the method of this invention has very small errors in reconstructing the differential phase, high accuracy, and fast algorithm speed, making it well-suited for multi-directional differential phase reconstruction in transverse shear interferometry.

[0253] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention without departing from its spirit and scope should be included within the protection scope of the present invention.

Claims

1. A multi-directional differential phase reconstruction method, characterized in that, The method includes the following steps: Step 1: Register the phase-shifted transverse shearing interferograms in the x and y directions respectively; Step 2: Obtain the differential phases in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm; Step 3: Transform the phase expansion into the optimal solution of the function, and use the least squares method of weighted iterative DCT to unwrap the differential phase and obtain the differential phase distribution in the x and y directions; Step 4: Perform edge detection and fitting on the phase-shifted transverse shearing interferograms in both the x and y directions, calculate the distance between the center points, and take the average value to obtain the shearing amount; specifically including: Step 4-1: Perform edge detection on the shearing interferogram, and then fill the gaps between the stripes; Step 4-2, perform edge detection again; Step 4-3: Fit two spot circles of the sheared beam from the left and right edges of the phase-shifted transverse shearing interference pattern in the x direction, respectively. Fit two spot circles of the sheared beam from the upper and lower edges of the phase-shifted transverse shearing interference pattern in the y direction, respectively. Calculate the distance between the centers of the two circles and take the average value to obtain the magnitude of the shearing. Steps 4-1 and 4-2 specifically utilize the Canny operator for edge detection; Step 5: Based on the differential phase distributions in the x and y directions obtained in Step 3 and the shearing amount obtained in Step 4, wavefront reconstruction is performed using the least squares method based on the difference Zernike polynomial; the specific process includes: According to the Zernike polynomial, the wavefront is represented as W(x,y): W(x,y)=∑a i Z i (x,y) Among them, Z i (x,y) represents the i-th Zernike polynomial, a i This indicates its corresponding coefficient; The shear wavefronts along the x and y directions are represented by ΔW, respectively. x (x,y) and ΔW y (x,y): ΔW x (x,y)=W(x+s,y)-W(x,y)=∑a xi [Z i (x+s,y)-Z i (x, y)]=∑a xi Z i (x,y) ΔW y (x,y)=W(x,y+s)-W(x,y)=∑a yi [Z i (x,y+s)-Z i (x, y)]=∑a yi Z i (x, y) Where s is the shear rate; The above relationship can be expressed in vector form: Right now In the formula, a x =[a x1 ,a x2 ,…,a xi ] T a y =[a y1 ,a y2 ,…,a yi ] T a x1 and a y1 T represents the Zernike coefficient of the i-th shear wavefront, respectively. x and T y The coefficient transformation matrix is ​​obtained by combining the relationship between the difference Zernike polynomials and the Zernike polynomials and eliminating the correlation columns. The Zernike coefficients before the measured wavefront are expressed as: Among them, T + It is the generalized inverse matrix of T; The wavefront under test is fitted with the Zernike coefficients of the wavefront under test, thus completing the reconstruction from the shear wavefront to the wavefront under test.

2. The multi-directional differential phase reconstruction method according to claim 1, characterized in that, Step 1, which involves registering the phase-shifted transverse shearing interferograms in the x and y directions, specifically includes: Step 1-1: Obtain the translation amount of each phase-shifted transverse shearing interferogram relative to the first phase-shifted transverse shearing interferogram using a phase-correlation image registration algorithm. Steps 1-2: Based on the translation amount, perform reverse translation on the phase-shifting transverse shearing interferograms to be matched so that the positions of each phase-shifting transverse shearing interferogram are matched.

3. The multi-directional differential phase reconstruction method according to claim 2, characterized in that, Step 1-1, the phase-correlation-based image registration algorithm, obtains the translation amount of each phase-shifted transverse shearing interferogram relative to the first phase-shifted transverse shearing interferogram, specifically including: Let f1 be the reference image of the phase-shifted transverse shearing interferogram, and f2 be the image to be registered that has been translated in the spatial domain. Then: f1 = f(x,y) f2 = f(x - Δx, y - Δy) In the formula, (Δx,Δy) represents the translation of the image to be registered f2 relative to the reference image f1 along the x-axis and y-axis, and f(x,y) represents the pixel value of the image at (x,y). Let the size of each phase-shifted transverse shearing interferogram be M×N. Performing Fourier transforms on f1 and f2 respectively, we have: F1(u,v)=FT{f(x,y)} In the formula, F1(u,v) and F2(u,v) are the spectrum diagrams after the Fourier transform of f1 and f2, respectively; Calculate the normalized power spectra of the two in the frequency domain and perform an inverse Fourier transform to obtain the unit impulse function C(x,y): In the formula, * represents the conjugate operation, and δ is the Dirac function; The translation amount can be obtained by determining the coordinates corresponding to the maximum value of the unit impulse function:

4. The multi-directional differential phase reconstruction method according to claim 3, characterized in that, Step 2 describes obtaining the differential phases enclosed in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm. Specifically, this includes: I1=A+Bcos(φ+π / 2)=A-Bsin(φ) I2=A+Bcos(φ+π)=A-Bcos(φ) l3=A+Bcos(φ+3π / 2)=A+Bsin(φ) I4=A+Bcos(φ+2π)=A+Bcos(φ) In the formula, I1, I2, I3, and I4 are the light intensities of the interference field, A is the background light intensity of the interference field, and B is the modulation index of the interference field.

5. The multi-directional differential phase reconstruction method according to claim 4, characterized in that, Step 3 involves transforming the phase expansion into solving for the optimal solution of the function. The least squares method of weighted iterative DCT is used to unwrap the differential phase, obtaining the differential phase distributions in the x and y directions. The specific process includes: Step 3-1: Represent the least squares problem as a vector form determined by multiple factors: Ax = b Further expressed using the cosine transform: A T Ax=A T b Where x is a phase value result vector of length MN, b is a wrapped phase difference vector of length N(M-1)+M(N-1), and T represents matrix transpose; Step 3-2, the weighted least squares problem is expressed as: WAx = Wb Further expressed using the discrete cosine transform, it can be represented as: A T W T WAx=A T W T Wb Where W is the weight matrix associated with the pixel weights; Order Q=A T W T WA, Then we have: Qφ=c in, is the weighted phase difference vector, and c is the weighted phase difference discrete phase Laplace vector; Decompose matrix Q into matrix P and difference matrix D; Step 3-3: Calculate c using the initial weight data and the wrapped phase difference vector; Steps 3-4: Set the initial conditions as iteration number k = 0, differential phase φ k =0, maximum number of iterations k max ; Steps 3-5: Iteratively calculate the differential phase φ k+1 : f k+1 =c-Dφ k Steps 3-6: Solve for ρ using the DCT least squares method. k : Rφ k+1 =ρ k In the formula, P is a vector φ k+1 The matrix to which the discrete Laplace operation is performed, ρ k It is a vector that contains a discrete Laplace operation on the phase difference of the package; Steps 3-7, determine k <k max If the condition is true, continue iterating from step 3-5 to step 3-6 until k = k. max Finally, we can solve for φ. k+1 .

6. A multi-directional differential phase reconstruction system based on the method of any one of claims 1 to 5, characterized in that, The system includes: The first module is used to perform position registration of the phase-shifted transverse shearing interferograms in the x and y directions respectively; The second module is used to obtain the differential phases wrapped in the x and y directions from the light intensity of the phase-shifted transverse shearing interferogram using a four-step phase-shifting algorithm. The third module is used to transform the phase expansion into the optimal solution of the solution function. Based on the least squares method of weighted iterative DCT, the differential phase is unwrapped to obtain the differential phase distribution in the x and y directions. The fourth module is used to perform edge detection and fitting on the phase-shifted transverse shearing interferograms in the x and y directions, calculate the distance between the centers of the circles and take the average value to obtain the shearing amount; The fifth module is used to perform wavefront reconstruction based on the differential phase distribution and shearing amount in the x and y directions obtained above, using the least squares method based on the differential Zernike polynomial.

7. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method according to any one of claims 1 to 5.

Images