A hybrid error correction phase recovery algorithm of frequency-shifted phase-shifted least square iteration

Abstract

The present application relates to a kind of variable frequency phase shift least square iteration hybrid error correction phase recovery algorithm.Based on least square iteration algorithm, by establishing a stripe model including time sequence light intensity fluctuation, phase shift error and Gamma distortion three error influences, introduce variable frequency stripe to provide minimum and sufficient number of numerical solution equation, for calculating the numerous to be solved parameters in the constructed stripe model.For the phase jump problem in some sampling point position in the solving process, due to the non-full rank of solving matrix, it is proposed to use regularization and combine local space Gamma distortion coefficient, background intensity and first-order harmonic intensity continuous and stripe intensity non-negative constraint to identify processing.The present application can realize the effective correction of mixed error including time sequence light intensity fluctuation, phase shift error and Gamma distortion in the case of using as few stripe diagram as possible, significantly improve phase demodulation precision.

Description

Technical Field

[0001] This invention relates to a high-precision phase recovery technique for phase-shifted fringe patterns, and in particular to a hybrid error correction phase recovery algorithm based on frequency-shifted least squares iteration, which belongs to the field of advanced optical detection technology. Background Technology

[0002] Phase-shifting algorithms offer advantages such as full-field measurement, high precision, point-by-point calculation, and good flexibility, making them a powerful fringe analysis tool in fields such as interferometry, holographic detection, fringe projection profilometry, and phase deflection. Traditional equal-interval phase-shifting algorithms rely on the captured phase-shifted fringe pattern strictly following a sine or cosine function distribution and require ideally equal step sizes for phase shift between frames. However, in actual measurements, various error factors often render these assumptions invalid, leading to phase demodulation errors, with wavy artifact errors being the most typical.

[0003] Among numerous influencing factors, the differences in fringe pattern intensity caused by light fluctuations in the light source or background, the phase shift error of the phase shifter, and the nonlinear Gamma distortion of optoelectronic devices are the main error sources that lead to reduced accuracy of wrapper phase demodulation during actual phase shift measurements. Therefore, how to simultaneously correct or eliminate the combined effects of the above-mentioned major error factors is one of the hot topics in the field of phase shift algorithm research, and it has significant research significance and practical application value.

[0004] To achieve more accurate phase retrieval, numerous scholars have proposed a series of phase-shifting algorithms. These include algorithms that only consider the influence of phase-shifting error factors, such as the classic Advanced Iterative Algorithm (AIA) and Principal Component Analysis (PCA) algorithms; and algorithms that only consider the influence of higher harmonic factors, such as the frequency-shifting least-squares iterative algorithm based on higher harmonic models and the fast combined frequency phase extraction algorithm. However, these algorithms only consider one of many error sources. When multiple error sources exist, the accuracy of the algorithm is severely limited, thus lacking versatility. While some scholars have conducted combination analysis and correction of pairwise errors based on the algorithms considering single error factors, existing phase-shifting algorithms still fail to comprehensively consider the combined influence of these three common error sources and effectively correct errors with as few fringe patterns as possible. This is of great significance for the research, application, and promotion of high-precision, high-performance, and universal phase-shifting algorithms. Summary of the Invention

[0005] This invention addresses the shortcomings of existing technologies by providing a universal phase-shifting algorithm that can recover the absolute phase of fringes with high precision using as few phase-shifting steps as possible. This algorithm effectively corrects mixed errors, including temporal intensity fluctuations, phase-shifting errors, and Gamma distortion, and significantly improves phase demodulation accuracy.

[0006] The technical solution to achieve the objective of this invention is to provide a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration, comprising the following steps:

[0007] Step 1: Initial value estimation: 1) The frequency conversion phase shift fringe pattern to be processed A phase-shifting algorithm based on different phase-shifting steps is used to process the data, obtaining initial values ​​of the parameters to be solved, including the error: wrapping phase. Background intensity First harmonic intensity coefficient Phase shift error Time-series light intensity fluctuation coefficient 2) The Gamma distortion coefficient γ in the fringe model (1) to be constructed 0 The initial value of (x,y) is set to 1; 3) Employ a phase unrolling algorithm to process the wrapped phase. Obtaining absolute phase

[0008] Step 2: Mixed error correction processing: 1) Construct a fringe model (1) including phase shift error, temporal intensity fluctuation and Gamma distortion error using the initial values ​​of the parameters to be solved obtained in step one: Where k is the iteration number; n = 0, 1, ..., N-1, N is the total number of phase shift steps; i = 0, 1, ..., L-1, L is the total number of temporal intensity fluctuation coefficients; m = 1, 2, ..., M, M represents the number of frequencies of the frequency conversion stripes; (x, y) represents the coordinates of the sampling point on the captured image plane; a mi Let A0(x,y) be the intensity fluctuation coefficient of the i-th time series of the fringe pattern at the m-th frequency; A0(x,y) is the background intensity; and B1(x,y) is the first harmonic intensity coefficient. The continuous phase distribution corresponding to the highest frequency fringe pattern; β m f m The ratio coefficient of the continuous phase distribution of the fringe pattern at a lower frequency to the continuous phase distribution of the fringe pattern at the highest frequency. f m ε is the frequency value of the fringe pattern at the m-th frequency; mn f mThe frequency and the phase shift amount of the fringe pattern at the nth step phase shift are determined and set as equal-step phase shifts; d mn f m Frequency, phase shift error of the fringe pattern at the nth phase shift; γ(x,y) is the Gamma distortion coefficient; 2) Divide the parameters to be solved in equation (1) into three groups, namely: Continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient γ(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient a mi The least squares method combined with regularization is used to perform step-by-step iterative calculations on the three sets of parameters to be solved. Step 1: Calculation of continuous phase distribution, background intensity, first harmonic intensity, and Gamma distortion coefficient: Step 1a: Construct matrix Y1 based on the difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model obtained under the k-th iteration, and determine the parameters to be solved A0, B1, γ and γ based on the fringe model equation (1). The gradient is used to construct matrix K1, and the least squares method combined with regularization is used to solve the matrix equation system (2): X1=(K1 T K1+ηE) -1 K1 T Y1 (2) Where I is the identity matrix and η is the regularization coefficient; In equation (3), the matrix X1 to be solved is A0, B1, γ and γ at the (k+1)th iteration and the kth iteration. The difference between them: In equation (4), matrix Y1 is the residual vector of the matrix equation system under the current step: In equation (5), matrix K1 represents the residual vectors for the parameters A0, B1, γ, and γ in equation (1) of the stripe model. The Jacobian matrix obtained by differentiation: in, For the stripe model (1), calculate the gradient with respect to the background intensity A0. To find the gradient of the fringe model equation (1) with respect to the first harmonic intensity B1, To calculate the gradient of the fringe model equation (1) with respect to the Gamma distortion coefficient γ, For the fringe model (1) with respect to the phase distribution Find the gradient; Then, the parameters to be solved in equation (6) are obtained from matrix X1 in equation (3). and Step 1b, the result calculated in Step 1a and The parameters are used to determine the transition points and to determine the background intensity of the f1 frequency stripe pattern. and first harmonic intensity The periodic threshold is used as the judgment condition: Where T n These are empirical values, and are generally set to be less than the peak intensity of the fringe pattern in the absence of noise. The value gradually decreases as the noise increases. When the judgment condition of equation (7) is met, the sampling point is determined to be an error residual point, and then the neighborhood is assigned according to equation (8) for the error residual point: Where u = 0, ±1, ..., ±U, U represents half of the range of values ​​of the neighboring points in the x-direction; v = 0, ±1, ..., ±V, V represents half of the range of values ​​of the neighboring points in the y-direction; Value(·) is the neighborhood assignment operation, which can be either the mean or interpolation. The parameters to be solved at the error residual point are obtained from Step 1b. as well as Step 2: Phase shift error calculation The result calculated in Step 1, which does not contain residual error, is the result of the k-th iteration. as well as The value of is used as the known value of subsequent Step 2 and Step 3 in the current iteration cycle to construct the fringe model (1) under the k-th iteration; the matrix Y2 is constructed based on the difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model (1) obtained under the k-th iteration, and the parameter d to be solved is determined based on the fringe model (1). mn The gradient is used to construct matrix K2, and the least squares method combined with regularization is used to solve the matrix equation system (9): X2=(K2 T K2+ηE) -1 K2 T Y2(9) In equation (10), the matrix X2 to be solved is the matrix d at the (k+1)th iteration and the kth iteration. mn The difference between them: In equation (11), matrix Y2 is the residual vector of the matrix equation system under the current step: In equation (12), matrix K2 represents the residual vector and the parameter d to be solved in equation (1) of the stripe model. mn The Jacobian matrix obtained by differentiation: in, The fringe model (1) represents the phase shift error d. mn Find the gradient; Then, the parameters to be solved in equation (13) are obtained from matrix X2 in equation (10). Step 3: Calculation of temporal light intensity fluctuation coefficient: The k-th iteration calculated by Step 2 The value of is used as the known value of the subsequent Step 3 in the current iteration cycle to construct the fringe model (1) under the k-th iteration; the matrix Y3 is constructed according to the difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model obtained under the k-th iteration, and the a is calculated according to the fringe model (1). mi The gradient of the parameters is used to construct matrix K3, and the least squares method combined with regularization is used to solve the matrix equation system of equation (14): X3=(K4 T K4+ηE) -1 K4 T Y4(14) In equation (15), the matrix X3 to be solved is the matrix a in the (k+1)th iteration and the matrix a in the kth iteration. mi The difference between them: In equation (16), matrix Y3 is the residual vector of the matrix equation system under the current step: In equation (17), matrix K3 represents the residual vector pair of the parameter a to be solved in equation (1) of the stripe model. mi The Jacobian matrix obtained by differentiation: in, For the fringe model equation (1), the temporal intensity fluctuation coefficient a mi Find the gradient; The parameters to be solved in equation (18) are obtained from matrix X3 in equation (15). After obtaining all unknown parameters through Steps 1, 2, and 3 above, a threshold judgment is performed to determine whether the iteration should terminate; the threshold calculation process is as follows (19): When the threshold condition is met or the number of iterations reaches the maximum set value, the iteration ends and the phase calculation result is output.

[0009] This invention provides a hybrid error correction phase recovery algorithm based on frequency-shifted phase-shift least squares iteration. In step one, the phase shift algorithm based on different phase shift steps is one of the following: an improved phase shift algorithm based on Lissajous ellipse fitting of two-step phase shift fringe patterns, or an improved phase shift algorithm based on advanced iterative algorithms of three-step or higher phase shift fringe patterns. In step one, the phase expansion method is a multi-frequency phase expansion method. In step two, the regularization process is one of Tikhonov regularization, total variational regularization, two-parameter shaping regularization, hybrid two-parameter regularization, and sparse structure constraint regularization, based on constraints of local spatial Gamma distortion coefficient, continuity of background intensity and first-order harmonic intensity, and non-negativity of fringe intensity.

[0010] This invention provides a hybrid error correction phase retrieval algorithm based on frequency-variable phase shift least squares iteration. The principle behind this algorithm is the introduction of frequency-variable fringes to address the inability of traditional single-frequency three-step or even two-step phase shift fringe algorithms to correct gamma distortion, and to obtain absolute phase information. Therefore, by using as few fringe patterns as possible, effective correction of hybrid errors, including phase shift error, temporal intensity fluctuations, and gamma distortion, can be achieved. Theoretically, at least a single-frequency four-step or dual-frequency two-step phase shift image can be used to correct hybrid errors, including phase shift error, temporal intensity fluctuations, and gamma distortion.

[0011] This invention provides a hybrid error correction phase recovery algorithm based on frequency-shift phase shift least squares iteration. It employs a method of constructing a fringe model that comprehensively considers the effects of three errors: phase shift error, temporal intensity fluctuation, and Gamma distortion, to approximate the frequency-shift fringe pattern to be processed. Frequency-shift fringes are introduced to provide a minimum and sufficient number of numerical solution equations for calculating the numerous parameters to be solved in the constructed fringe model. Utilizing the fixed ratio between the absolute phases of each frequency-shift fringe and the property that the degree of Gamma distortion is consistent under the same measurement system environment, the three error influencing factors are grouped and iteratively calculated in steps for correction. For the phase jump problem in some regions caused by the solution matrix not being full rank at certain sampling points during the solution process, a regularization method is proposed, combined with constraints such as the continuity of local spatial Gamma distortion coefficients, background intensity, and first-order harmonic intensity, as well as the non-negativity of fringe intensity. This improves the numerical stability when solving for non-full-rank matrices, thereby demodulating high-precision absolute phase distribution results.

[0012] Compared with the prior art, the significant advantages of the present invention are: 1. Under the same test conditions, the least number of fringe patterns can be used to achieve effective correction of mixed errors, including temporal intensity fluctuations, phase shift errors and Gamma distortion, through least squares method and grouped distribution iterative calculation, significantly improving the phase demodulation accuracy; 2. A method is proposed to use regularization and combine it with constraints such as the continuity of local spatial Gamma distortion coefficient, background intensity and first harmonic intensity, and non-negativity of fringe intensity to identify phase jump points, thereby improving the numerical stability of the least squares iterative method for solving non-full rank matrices. Attached Figure Description

[0013] Figure 1 The flowchart illustrates a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration, provided in an embodiment of the present invention.

[0014] Figure 2 This is a schematic diagram illustrating the specific algorithm iteration process of the hybrid error correction processing in a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration provided in an embodiment of the present invention.

[0015] Figure 3The following are simulation results comparisons of a hybrid error correction phase recovery algorithm based on frequency-shifting least squares iteration provided in this embodiment of the invention. When processing a dual-frequency dual-step phase-shifting fringe pattern containing the effects of hybrid errors including temporal intensity fluctuations, phase shift errors, and Gamma distortion, (a) shows the object phase calculation result of the traditional dual-step phase-shifting algorithm LEF, (b) shows the object phase calculation result after initial value estimation, (c) shows the object phase calculation result after hybrid error correction, (d) shows the error residue of the phase calculation result of the traditional dual-step phase-shifting algorithm LEF, (e) shows the error residue of the phase calculation result after initial value estimation, and (f) shows the error residue of the phase calculation result after hybrid error correction.

[0016] Figure 4 A comparison of phase error cross-sections from simulation results of a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration provided in an embodiment of the present invention. Detailed Implementation

[0017] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0018] Example 1 See appendix Figure 1 This is a flowchart of a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration provided in this embodiment. Subsequent embodiments will use a dual-frequency, dual-step phase shift fringe pattern as an example. The actual acquired dual-frequency, dual-step phase shift fringe pattern will be processed in the following two steps:

[0019] Step 1: Initial value estimation: (1) The Lissajous ellipse fitting (LEF) optimization iterative algorithm was used to process the frequency conversion phase shift fringe patterns to be processed, and the initial values ​​of the parameters to be solved with errors were calculated: wrapping phase Background intensity First harmonic intensity coefficient Phase shift error and temporal light intensity fluctuation coefficient The processing steps of the LEF optimization iterative algorithm are as follows: First, a fringe model is established that only considers phase shift error and temporal intensity fluctuation: Where, k' is the iteration number; n = 0, 1, ..., N-1, N = 2, is the total number of phase shift steps; i = 0, 1, ..., L-1, L = 2, is the total number of temporal intensity fluctuation coefficients; m = 1, 2, ..., M, M = 2, represents the number of frequencies of the frequency conversion stripes; (x, y) represents the coordinates of the sampling point on the captured image plane; a' miφ' is the intensity fluctuation coefficient of the i-th time series light intensity in the fringe pattern at the m-th frequency; A'0(x,y) is the background intensity; B'1(x,y) is the first harmonic intensity coefficient; φ' m (x,y) represents the continuous phase distribution corresponding to the fringe pattern; ε mn f m The frequency and the phase shift amount used to determine the fringe pattern at the nth phase shift are generally set as equal-step phase shifts; d' mn Let be the phase shift error of the fringe pattern at the m-th frequency and the n-th phase shift. Initial values ​​for the parameters in the above fringe model are first estimated, as shown in Step 0 below.

[0020] Step 0 LEF initial value estimation: The temporal light intensity fluctuation coefficient was calculated using polynomial fitting. Where mean(·) is the mean operation. The image shows the dual-frequency, dual-step phase-shifting fringe pattern to be processed.

[0021] Perform the operations of formulas (3) and (4) on the actual acquired phase shift fringe pattern: The above obtained and Ellipse fitting was performed on the two sets of data: The LEF algorithm was used to find the center coordinates x of the ellipse in formula (5). m0 y m0 and major and minor axes a mx a my The parameters are calculated, and then the wrapping phase is calculated based on these four parameters. The temporal intensity fluctuation coefficient is obtained through formulas (2) and (6). and package phase In the case of background intensity, the following formulas (7) to (11) are used to calculate the background intensity respectively. First harmonic intensity coefficient and phase shift error Initial value: Where N x N represents the maximum number of sampling points along the x-axis on the captured image plane. yThis represents the maximum number of sampling points along the y-axis on the captured image plane.

[0022] A'0(x,y)=mean(C1)·ones(N x N y (9) Among them ones(N) x N y ) represents constructing an N x Line N y A matrix of all identical columns.

[0023] The above calculation yields the package phase. Background intensity First harmonic intensity coefficient Phase shift error and temporal light intensity fluctuation coefficient Substitute into formula (1) to construct the stripe model. Then, based on the stripe model of formula (1), carry out the subsequent iterative calculations in Steps 1, 2, and 3.

[0024] Step 1: Calculation of continuous phase distribution, background intensity, and first harmonic intensity: In the phase shift fringe pattern at the m-th frequency, the parameters to be solved are A'0, B'1, and... It has the property of being approximately equal. The solution formula is shown in equation (12). The least squares method is used in combination with regularization. This linear system of equations has two equations and involves two unknowns.

[0025] X'1=(K'1 T K'1+ηE) -1 K'1 T Y'1 (12) Where E is the identity matrix and η is the regularization coefficient, which is an empirical parameter.

[0026] Matrix X'1 represents the difference between A'0 and B'1 at the (k'+1)th iteration and the k'th iteration, and is the matrix to be solved: Matrix Y'1 is the residual vector of the solution matrix in this step: Matrix K'1 represents the Jacobian matrix generated by differentiating the residual vector with respect to the parameters A'0 and B'1 in formula (1) of the stripe model formula: in For the stripe model equation (1), calculate the gradient with respect to the background intensity A'0. Find the gradient of the first harmonic B'1 for the fringe model equation (1).

[0027] The parameters for the next iteration are obtained using matrix X'1: Obtain the parameters to be solved in the k'th iteration. Then, the results are obtained through equations (17-19).

[0028] Step 2: Phase shift error calculation In the phase shift fringe pattern at the m-th frequency, the parameter d' to be solved is... mn The values ​​between them are not equal. The solution formula is shown in equation (20). The least squares method is used in combination with regularization. This linear system of equations has two equations and involves two unknowns.

[0029] X'2=(K'2 T K'2+ηE) -1 K'2 T Y'2 (20) Matrix X'2 represents d' in the (k'+1)th iteration and the k'th iteration. mn The difference between them is the matrix to be solved: Matrix Y'2 is the residual vector of the solution matrix in this step: Matrix K'2 represents the residual vector pair, and the parameter d' to be solved in the formula for the stripe model is d'. mn The Jacobian matrix generated by differentiation: in For the fringe model (1), the phase shift error d' mn Find the gradient.

[0030] The parameters for the next iteration are obtained using matrix X'2:

[0031] Step 3: Calculation of temporal light intensity fluctuation coefficient: In the phase shift fringe pattern at the m-th frequency, the parameter a' to be solved is... mi The values ​​between them are not equal. The solution formula is shown in equation (25). The least squares method is used to solve the system of linear equations, which has two equations and involves two unknowns.

[0032] X'3=(K'3 T K'3+ηE) -1 K'3 T Y'3 (25) Matrix X'3 represents a' in the (k'+1)th iteration and the k'th iteration. mi The difference between them is the matrix to be solved: Matrix Y'3 is the residual vector of the solution matrix in this step: Matrix K'3 represents the residual vector pair, and the parameter a' to be solved in the formula for the stripe model is... mi The Jacobian matrix generated by differentiation: in For the fringe model equation (1), the temporal intensity fluctuation coefficient a' mi Find the gradient.

[0033] The parameters for the next iteration are obtained using matrix X'3: After obtaining all unknown parameters through the above three steps, a threshold check is performed to determine whether the iteration should terminate. When the threshold condition is met or the number of iterations reaches the maximum set value, the iteration ends and the phase calculation result is output. The threshold calculation process is shown below (thresholds Δ'1 and Δ'2 are empirically set to 5 × 10⁻⁶). -6 and 5×10 -6 ): After satisfying the above threshold conditions, the following initial value is obtained: Wrapping Phase Background intensity First harmonic intensity coefficient Phase shift error and temporal light intensity fluctuation coefficient

[0034] (2) Set the Gamma distortion coefficient γ in the following formula (31) 0 (x,y) is initially set to 1; (3) Using the phase unwrapping algorithm, wrap the phase Unfolding into absolute phase

[0035] Step 2: Mixed error correction processing: (1) Construct a fringe model that includes the effects of three types of errors: phase shift error, temporal intensity fluctuation and Gamma distortion, using the initial values ​​of the parameters to be solved obtained in step one; (2) Divide the parameters to be solved in the stripe model into three groups: Continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient γ(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient a mi The least squares method combined with regularization is used for step-by-step iterative calculation; the calculation ends after multiple iterations until the convergence condition is met, and the high-precision absolute phase and the above parameters are output.

[0036] See appendix Figure 2 This is a flowchart of the calculation process for hybrid error correction in a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration provided in this embodiment.

[0037] (1) Construct a fringe model with the initial values ​​of the parameters to be solved obtained in step one, including the effects of three types of errors: phase shift error, time-series light intensity fluctuation and higher harmonics, as shown in equation (29): Where k is the iteration number; n = 0, 1, ..., N-1, N = 2, is the total number of phase shift steps; i = 0, 1, ..., L-1, L = 2, is the total number of temporal intensity fluctuation coefficients; m = 1, 2, ..., M, M = 2, represents the number of frequencies of the frequency conversion stripes; (x, y) represents the coordinates of the sampling point on the captured image plane; a mi Let A0(x,y) be the intensity fluctuation coefficient of the i-th time series of the fringe pattern at the m-th frequency; A0(x,y) is the background intensity; and B1(x,y) is the first harmonic intensity coefficient. The continuous phase distribution corresponding to the highest frequency fringe pattern; β m f m The ratio coefficient of the continuous phase distribution of the fringe pattern at a lower frequency to the continuous phase distribution of the fringe pattern at the highest frequency. f m ε is the frequency value of the fringe pattern at the m-th frequency; mn f m The frequency and the phase shift amount used to determine the fringe pattern at the nth phase shift are generally set as equal-step phase shifts; d mn f m Frequency, phase shift error of the fringe pattern at the nth phase shift; γ(x,y) is the Gamma distortion coefficient; (2) Divide the parameters to be solved in equation (31) into three groups: Continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient γ(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient a mi The least squares method combined with regularization is used to perform step-by-step iterative calculations on the three sets of parameters to be solved, including the following steps: Step 1: Calculation of continuous phase distribution, background intensity, first harmonic intensity, and Gamma distortion coefficient: In a dual-frequency, dual-step phase-shift fringe pattern, the parameters to be solved are A0, B1, γ, and... They possess approximately equal characteristics. The difference between the frequency-shifted phase-shift fringe pattern to be processed and the fringe model under the k-th iteration, as well as the fringe model's relationship to A0, B1, γ, and... The gradient of the parameters is used to construct a set of solution matrix equations that continuously approximate the true values ​​of the parameters in each part of the dual-frequency dual-step phase-shifted fringe pattern to be processed. The least squares method is then used in conjunction with regularization processing, as shown in equation (32). This linear equation set has 4 equations and involves 4 unknowns.

[0038] X1=(K1 T K1+ηE) -1 K1 T Y1 (32) Where E is the identity matrix and η is the regularization coefficient, which is an empirical parameter.

[0039] Matrix X1 represents A0, B1, γ, and... in the (k+1)th and kth iterations. The difference between them is the matrix to be solved: Matrix Y1 is the residual vector of the solution matrix in this step: Matrix K1 represents the residual vector pair of the parameters to be solved in the stripe model formula, namely A0, B1, γ, and... The Jacobian matrix generated by differentiation: in To calculate the gradient of the stripe model equation (31) with respect to the background intensity A0, To calculate the gradient of the fringe model (31) with respect to the first harmonic B1, To calculate the gradient of the fringe model equation (31) with respect to the Gamma distortion coefficient γ, For the fringe model (31) phase distribution Find the gradient; The parameters for the next iteration are obtained using matrix X1: The parameters to be solved in the k-th iteration are obtained through the calculation in Step 1. as well as As known values ​​for the next two steps within the current iteration cycle.

[0040] Step 2: Phase shift error calculation In a dual-frequency, dual-step phase-shift fringe pattern, the parameter d to be solved is... mn The values ​​between them are not equal. The difference between the frequency-shifted phase-shift fringe pattern to be processed and the fringe model under the k-th iteration is used, along with the fringe model's influence on the d... mn The gradient of the parameters is used to construct a set of solution matrix equations that continuously approximate the true values ​​of the parameters in each part of the dual-frequency dual-step phase-shifted fringe pattern to be processed. The least squares method is then used in conjunction with regularization processing, as shown in equation (37). This linear equation set has 4 equations and involves 4 unknowns.

[0041] X2=(K2 T K2+ηE) -1 K2 T Y2(37) Matrix X2 represents d during the (k+1)th iteration and the kth iteration. mn The difference between them is the matrix to be solved: Matrix Y2 is the residual vector of the solution matrix in this step: Matrix K2 represents the residual vector pair, which is the parameter d to be solved in the formula for the stripe model. mn The Jacobian matrix generated by differentiation: in The fringe model (31) represents the phase shift error d. mn Find the gradient.

[0042] Then, the parameters for the next iteration are obtained using matrix X2: The parameters to be solved in the k-th iteration are obtained through the calculation in Step 2. As a known value for the next step within the current iteration cycle.

[0043] Step 3: Calculation of temporal light intensity fluctuation coefficient: In a dual-frequency, dual-step phase-shift fringe pattern, the parameter a to be solved is... mi The values ​​between them are not equal. The difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model under the k-th iteration is used, along with the fringe model's influence on the value of a. miThe gradient of the parameters is used to construct a set of solution matrix equations that continuously approximate the true values ​​of the parameters in each part of the dual-frequency, dual-step phase-shifted fringe pattern to be processed. The matrix is ​​then solved by combining regularization processing, as shown in equation (42). This linear equation set has 4 equations and involves 4 unknowns.

[0044] X3=(K3 T K3+ηE) -1 K3 T Y3 (42) Matrix X3 represents the values ​​of a in the (k+1)th iteration and the kth iteration. mi The difference between them is the matrix to be solved: Matrix Y3 is the residual vector of the solution matrix in this step: Matrix K3 represents the residual vector pair, which is the parameter to be solved in the formula for the stripe model. mi The Jacobian matrix generated by differentiation: in The fringe model (31) represents the temporal intensity fluctuation coefficient a. mi Find the gradient.

[0045] The parameters for the next iteration are obtained using matrix X3: The parameters to be solved at the k-th iteration number are obtained through the calculation in Step 3.

[0046] After obtaining all unknown parameters through the above three steps, a threshold judgment is performed to determine whether the iteration should terminate. When the threshold condition is met or the number of iterations reaches the maximum set value, the iteration ends and the phase calculation result is output. The threshold calculation process is shown below (thresholds Δ1, Δ2, and Δ3 are empirically set to 1×10⁻⁶). -6 1×10 -6 and 1×10 -6 ): Step a. Handling residual error points: Given the non-uniform spatial distribution of fringe intensity in the dual-frequency, dual-step fringe pattern to be processed, and the significant residual error in Step 1a, a mechanism for identifying residual error points is introduced to address the issues arising from the calculations in Step 1a. and γ k+1 The parameters are used to determine the transition point and adjust the background modulation of the f1 frequency stripe pattern. and first harmonic modulation The periodic threshold is used as the judgment condition: Where T n These are empirical values, determined based on the intensity of the actual stripe pattern. When the above judgment conditions are met, the sampling point is determined to be an error residual point, and an appropriate value will be assigned to the neighborhood of the error residual point: Where u = 0, ±1, ..., ±U, U represents half of the range of values ​​for the neighboring points in the x-direction; v = 0, ±1, ..., ±V, V represents half of the range of values ​​for the neighboring points in the y-direction; Value(·) is the value operation in the neighborhood, which can be the mean, median, or interpolation, etc. Through the above operations, residual error points are significantly corrected, achieving high-precision phase calculation.

[0047] See appendix Figure 3 The figures show a comparison of simulation results for a hybrid error correction phase recovery algorithm based on frequency-shifted phase-shift least squares iteration provided by an embodiment of the present invention. When processing a dual-frequency, dual-step phase-shifted fringe pattern containing the effects of hybrid errors including temporal intensity fluctuations, phase shift errors, and Gamma distortion, (a) shows the object phase calculation result of the traditional LEF algorithm, (b) shows the object phase calculation result of the improved LEF algorithm, (c) shows the object phase calculation result of the proposed algorithm; (d) shows the residual error of the phase calculation result of the traditional LEF algorithm (RMSE = 0.5465, Max error = 1.2192), (e) shows the residual error of the phase calculation result of the improved LEF algorithm (RMSE = 0.1650, Max error = 0.6997), and (f) shows the residual error of the phase calculation result of the proposed algorithm (RMSE = 0.0006, Max error = 0.0018). It is evident that the traditional LEF algorithm cannot address the impact of errors such as temporal intensity fluctuations, phase shift errors, and gamma distortion. Consequently, the calculated object phase distribution exhibits significant residual periodic fringe artifacts. The improved LEF algorithm, which only considers temporal intensity fluctuations and phase shift errors, also retains noticeable periodic fringe artifacts in the object phase distribution after processing fringe patterns that incorporate these factors. The proposed algorithm comprehensively considers the combined effects of temporal intensity fluctuations, phase shift errors, and gamma distortion, resulting in a high degree of consistency between the calculated object phase and the actual object phase, achieving high-precision phase calculation.

[0048] See appendix Figure 4This figure shows a comparison of the phase error cross-sections of the simulation results of a hybrid error correction phase recovery algorithm based on frequency-shift phase-shift least squares iteration, provided by an embodiment of the present invention. The figure compares the residual error cross-sections of the object phase calculation results of the traditional LEF algorithm, the improved LEF algorithm, and the proposed algorithm. Analysis shows that when processing dual-frequency, dual-step phase-shift fringe patterns containing the effects of mixed errors such as temporal intensity fluctuations, phase shift errors, and Gamma distortion, the phase calculation accuracy of the proposed algorithm is significantly better than that of the traditional and improved LEF algorithms.

[0049] Example 2 This embodiment provides a hybrid error correction phase recovery algorithm based on frequency-shifted phase-shift least squares iteration. Taking a dual-frequency, dual-step fringe pattern as an example, the method provided in Embodiment 1 is adopted, including the following steps: Step 1: Initial value estimation: (1) The improved LEF algorithm is used to process the dual-frequency dual-step phase-shift fringe patterns to be processed, and the initial values ​​of the parameters to be solved with errors are calculated: wrapping phase Background intensity First harmonic intensity coefficient Phase shift error and temporal light intensity fluctuation coefficient (2) Set the Gamma distortion coefficient The initial value is 1; (3) Using the phase unrolling algorithm, the wrapped phase is... Unfolding into absolute phase Step 2: Mixed error correction processing: (1) Construct a fringe model including the effects of three types of errors: phase shift error, temporal intensity fluctuation, and Gamma distortion, based on the initial values ​​of the parameters to be solved obtained in Step 1; (2) Divide the parameters to be solved in the fringe model into three groups: Continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient γ(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient a mi The least squares method combined with regularization is used for step-by-step iterative calculation; the calculation ends after multiple iterations until the convergence condition is met, and the high-precision absolute phase and the above parameters are output.

[0050] In the hybrid error correction phase recovery algorithm for frequency-shift phase shift least squares iteration provided in this embodiment, the hybrid error correction processing in step two uses the parameters to be solved in the frequency-shift phase shift fringe pattern calculated in step one as initial values ​​to construct a fringe model including the effects of three types of errors: phase shift error, temporal intensity fluctuation, and Gamma distortion; then, the parameters to be solved in the fringe model are divided into three groups: Continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient γ(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient a mi The least squares method is used for iterative calculation, and regularization is adopted in combination with the constraint that the local spatial Gamma distortion coefficient is approximately unchanged and the fringe intensity is non-negative to deal with the error residual points.

[0051] In the hybrid error correction phase recovery algorithm for frequency-shift phase shift least squares iteration provided in this embodiment, step one, based on the improved LEF algorithm, calculates the unknown parameters in each frequency fringe pattern after considering the influence of phase shift error and timing light intensity fluctuations; wherein the unknown parameters are divided into three groups for iterative calculation: Encapsulate the phase distribution φ1(x,y), background intensity A0(x,y), and first harmonic intensity coefficient B1(x,y); Phase shift error d mn ; Temporal intensity fluctuation coefficient α mi This embodiment provides a hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration. The phase shift algorithm based on different phase shift steps in sub-step (1) of step one is either an improved phase shift algorithm based on Lissajous ellipse fitting based on two-step phase shift fringe pattern or an improved phase shift algorithm based on advanced iterative algorithm based on three-step or more phase shift fringe pattern.

[0052] This embodiment provides a hybrid error correction phase recovery algorithm for frequency-shifted phase-shift least squares iteration, wherein the phase expansion technique in sub-step (3) of step one is the multi-frequency phase expansion method.

[0053] This embodiment provides a hybrid error correction phase recovery algorithm for frequency-shifted phase shift least squares iteration. The regularization process in sub-step (2) of step two is one of Tikhonov regularization, total variational regularization, two-parameter shaping regularization, hybrid two-parameter regularization and sparse structure constraint regularization, combined with constraints of local spatial Gamma distortion coefficient, continuous background intensity and first harmonic intensity, and non-negative fringe intensity.

Claims

1. A hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration, characterized in that... Includes the following steps: Step 1: Initial value estimation: 1) The frequency conversion phase shift fringe pattern to be processed A phase-shifting algorithm based on different phase-shifting steps is used to process the data, obtaining initial values ​​of the parameters to be solved, including the error: wrapping phase. Background intensity 1st harmonic intensity coefficient Phase shift error Time-series light intensity fluctuation coefficient ; 2) The Gamma distortion coefficient in the stripe model (1) to be constructed The initial value is set to 1; 3) Employ a phase unrolling algorithm to process the wrapped phase. Obtaining absolute phase ; Step 2: Mixed error correction processing: 1) Construct a fringe model (1) with the initial values ​​of the parameters to be solved obtained in step one, including phase shift error, temporal intensity fluctuation and Gamma distortion error: (1) Where k is the iteration number; n = 0, 1, ..., N-1, N is the total number of phase shift steps; i = 0, 1, ..., L-1, L is the total number of temporal intensity fluctuation coefficients; m = 1, 2, ..., M, M represents the number of frequencies of the frequency conversion stripes; (x, y) represents the coordinates of the sampling point on the captured image plane; a mi Let A0(x,y) be the intensity fluctuation coefficient of the i-th time series of the fringe pattern at the m-th frequency; A0(x,y) is the background intensity; and B1(x,y) is the first harmonic intensity coefficient. The continuous phase distribution corresponding to the highest frequency fringe pattern; f m The ratio coefficient of the continuous phase distribution of the fringe pattern at a lower frequency to the continuous phase distribution of the fringe pattern at the highest frequency. f m The frequency value of the fringe pattern at the m-th frequency; f m Frequency, the phase shift amount of the fringe pattern at the nth phase shift, set as a constant step phase shift; d mn f m Frequency, phase shift error of the fringe pattern at the nth phase shift; Gamma distortion coefficient; 2) Divide the parameters to be solved in equation (1) into three groups, namely: a) continuous phase distribution Background intensity A0(x,y), first harmonic intensity coefficient B1(x,y), and Gamma distortion coefficient b. Phase shift error d. mn c. Temporal intensity fluctuation coefficient a mi The least squares method combined with regularization is used to perform step-by-step iterative calculations on the three sets of parameters to be solved. Step 1: Calculation of continuous phase distribution, background intensity, first harmonic intensity, and Gamma distortion coefficient: Step 1a: Construct matrix Y1 based on the difference between the frequency-shifted fringe pattern to be processed and the fringe model obtained under the kth iteration, and construct matrix K1 based on the gradients of the parameters A0, B1, γ and φ1 to be solved according to the fringe model equation (1). Solve the matrix equation system of equation (2) using the least squares method combined with regularization: (2) Where E is the identity matrix and η is the regularization coefficient; In equation (3), the matrix X1 to be solved is the difference between A0, B1, γ and φ1 in the (k+1)th iteration and the kth iteration: (3) In equation (4), matrix Y1 is the residual vector of the matrix equation system under the current step: (4) Equation (5) matrix K1 is the Jacobian matrix obtained by differentiating the residual vector with respect to the unsolved parameters A0, B1, γ and φ1 in the stripe model equation (1): (5) in, For the stripe model (1), calculate the gradient with respect to the background intensity A0. To calculate the gradient of the fringe model equation (1) with respect to the first harmonic intensity B1, For the stripe model equation (1), calculate the gradient of the Gamma distortion coefficient γ. For the fringe model (1), find the gradient of the phase distribution φ1; Then, the parameters to be solved in equation (6) are obtained from matrix X1 in equation (3). , , and : (6) Step 1b, the result calculated in Step 1a , , and The parameters are used to determine the transition points, and the background intensity of the f1 frequency value fringe pattern is determined. and first harmonic intensity The periodic threshold is used as the judgment condition: (7) Where T n This is an empirical value, set to be smaller than the peak intensity of the fringe pattern in the absence of noise, and the value gradually decreases as the noise increases. When the judgment condition of equation (7) is met, the sampling point is determined to be an error residual point, and then the neighborhood is assigned according to equation (8) for the error residual point: (8) Where u = 0, ±1, ..., ±U, U represents half of the range of values ​​of the neighboring points in the x-direction; v = 0, ±1, ..., ±V, V represents half of the range of values ​​of the neighboring points in the y-direction; Value(·) is the neighborhood assignment operation, which can be either the mean or interpolation. The parameters to be solved at the error residual point are obtained from Step 1b. , , as well as ; Step 2: Phase shift error calculation: The result calculated in Step 1, which does not contain residual error, is the result of the k-th iteration. , , as well as The value of is used as the known value of subsequent Step 2 and Step 3 in the current iteration cycle to construct the fringe model (1) under the k-th iteration; the matrix Y2 is constructed based on the difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model (1) obtained under the k-th iteration, and the parameter d to be solved is determined based on the fringe model (1). mn The gradient is used to construct matrix K2, and the least squares method combined with regularization is used to solve the matrix equation system (9): (9) In equation (10), the matrix X2 to be solved is the matrix d at the (k+1)th iteration and the kth iteration. mn The difference between them: (10) In equation (11), matrix Y2 is the residual vector of the matrix equation system under the current step: (11) In equation (12), matrix K2 represents the residual vector and the parameter d to be solved in equation (1) of the stripe model. mn The Jacobian matrix obtained by differentiation: (12) in, The fringe model (1) represents the phase shift error d. mn Find the gradient; Then, the parameters to be solved in equation (13) are obtained from matrix X2 in equation (10). : (13) Step 3: Calculation of temporal light intensity fluctuation coefficient: The k-th iteration calculated by Step 2 The value of is used as the known value of the subsequent Step 3 in the current iteration cycle to construct the fringe model (1) under the k-th iteration; the matrix Y3 is constructed according to the difference between the frequency conversion phase shift fringe pattern to be processed and the fringe model obtained under the k-th iteration, and the a is calculated according to the fringe model (1). mi The gradient of the parameters is used to construct matrix K3, and the least squares method combined with regularization is used to solve the matrix equation system (14): (14) In equation (15), the matrix X3 to be solved is the matrix a at the (k+1)th iteration and the kth iteration. mi The difference between them: (15) In equation (16), matrix Y3 is the residual vector of the matrix equation system under the current step: (16) In equation (17), matrix K3 represents the residual vector pair of the parameter a to be solved in equation (1) of the stripe model. mi The Jacobian matrix obtained by differentiation: (17) in, For the fringe model (1), the temporal intensity fluctuation coefficient a mi Find the gradient; The parameters to be solved in equation (18) are obtained from matrix X3 in equation (15). : (18) After obtaining all unknown parameters through Steps 1, 2, and 3 above, a threshold judgment is performed to determine whether the iteration should terminate; the threshold calculation process is as follows (19): (19) When the threshold condition is met or the number of iterations reaches the maximum set value, the iteration ends and the phase calculation result is output. .

2. The hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration according to claim 1, characterized in that: In step one, the phase shift algorithm based on different phase shift steps is one of the improved phase shift algorithm based on Lissajous ellipse fitting based on two-step phase shift fringe patterns and the improved phase shift algorithm based on advanced iterative algorithms based on three-step or higher phase shift fringe patterns.

3. The hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration according to claim 1, characterized in that: In step one, the phase expansion algorithm is a multi-frequency phase expansion method.

4. The hybrid error correction phase recovery algorithm based on frequency conversion phase shift least squares iteration according to claim 1, characterized in that: In step two, the regularization process is one of Tikhonov regularization, total variational regularization, two-parameter shaping regularization, hybrid two-parameter regularization, and sparse structure constraint regularization, based on the constraints of local spatial Gamma distortion coefficient, continuous background intensity and first harmonic intensity, and non-negative fringe intensity.

Images