Three-dimensional measurement method for non-equal period phase shift

Abstract

The invention discloses a non-equal period phase shift phase unwrapping method, which comprises the following steps of: 1, constructing a wavelength-De Bruijn element mapping table, and generating a non-equal period phase shift stripe group; 2, projecting the phase shift fringes to a measured object frame by frame, obtaining a modulation pattern, and decoding the modulation pattern to obtain a non-linear wrapped phase; 3, constructing an approximate phase to carry out local phase unwrapping on the nonlinear wrapped phase to obtain a nonlinear local phase; and 4, calculating a phase order by minimizing a nonlinear error on the nonlinear local phase, performing phase unwrapping in combination with a nonlinear wrapped phase, and obtaining a linear absolute phase of the measured object through piecewise linear mapping. According to the measurement method, high-precision absolute phase calculation can be realized only by using three projection patterns, and the three-dimensional measurement efficiency is greatly improved.

Description

Technical Field

[0001] The invention relates to the technical field of three-dimensional shape measurement, and in particular to a phase unwrapping method of non-equiperiodic phase shift. Background Art

[0002] With the advancement of image processing, optics and computer technology, optical image processing technology based on optoelectronic imaging theory has developed rapidly. In the past few decades, this technology has been widely studied and applied, and has also led to major breakthroughs in computer vision technology. As people's demand for information richness and diversity continues to increase, the evolution of computer vision technology has gradually expanded from two-dimensional images to three-dimensional space. Therefore, how to efficiently obtain high-quality spatial stereoscopic visual perception information has become one of the current research focuses in the field of computer vision.

[0003] As a low-cost, high-precision and efficient non-contact measurement method, optical 3D measurement technology is widely used in robot navigation, industrial modeling, microscopic measurement and other fields. The mainstream structured light measurement method includes fringe projection profile measurement, which uses a projector to project a grating pattern carrying phase information, and a camera synchronously captures the image modulated by the object being measured, and uses the phase shift method to decode it to obtain the phase information. Since the contour information of the object surface is hidden in the phase information, and the phase information can only be obtained through the arc tangent function, the effective range of the phase is limited to (-π,π], resulting in phase ambiguity. To eliminate the phase ambiguity, phase unwrapping technology is used to restore the continuous absolute phase.

[0004] To robustly unfold the wrapped phase into the absolute phase, it is necessary to use phase unwrapping technology to perform phase unwrapping. After decades of development, researchers have proposed many effective phase unwrapping methods, which can be divided into spatial phase unwrapping and temporal phase unwrapping. The spatial phase unwrapping method does not require the use of additional modes. It usually uses the phase relationship of spatial neighborhood pixels to achieve phase unwrapping. However, when the field of view contains depth mutations or isolated objects, phase unwrapping often makes mistakes. The temporal phase unwrapping method can perfectly solve the limitations of the spatial phase unwrapping method. It projects additional encoding information on the scene to uniquely identify each phase order of the wrapped phase, so there is no error propagation in the spatial domain. However, the additional multi-frame data acquisition of the temporal phase unwrapping method not only limits the measurement speed, but also accumulates time domain errors and produces additional artifacts. Therefore, an important goal is a reliable phase unwrapping method that does not require additional modes. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a phase unwrapping method for non-uniform periodic phase shift.

[0006] The technical solution of the invention is as follows:

[0007] A three-dimensional measurement method of non-uniform periodic phase shift, comprising the following steps:

[0008] Step 1: Construct a wavelength-De Bruijn element mapping table to generate a non-uniform periodic phase-shifted fringe group;

[0009] Step 2: Project the phase-shifted fringes onto the object under test frame by frame, obtain the modulation pattern and decode it to obtain the nonlinear wrapped phase;

[0010] Step 3: Construct an approximate phase and perform a local phase unwrapping operation on the nonlinear wrapped phase to obtain a nonlinear local phase;

[0011] Step 4: Calculate the phase order by minimizing the nonlinear error on the nonlinear local phase, then combine the phase order with the nonlinear wrapped phase for phase unfolding, and then obtain the linear absolute phase of the object under test through piecewise linear mapping to complete the three-dimensional measurement.

[0012] Step 1 includes:

[0013] Step 1.1, select m wavelengths and De Bruijn elements, and construct a wavelength-De Bruijn element mapping table;

[0014] The wavelength-De Bruijn element mapping table is:

[0015] De Bruijn Element <![CDATA[db1]]> <![CDATA[db2]]> ... <![CDATA[db m ]]> wavelength <![CDATA[f1]]> <![CDATA[f2]]> ... <![CDATA[f m ]]>

[0016] Step 1.2, generating an m-element n-order De Bruijn wavelength sequence and generating a non-uniform periodic phase shift code;

[0017] The De Bruijn wavelength sequence is:

[0018] F DB ={f db_1 ,f db_2 ,f db_3 ,...f db_mn}

[0019] The non-equal period phase shift encoding formula in step 1.2 is:

[0020]

[0021] in

[0022]

[0023] is the phase shift pattern, x p and pdenote the horizontal and vertical coordinates of the projector plane, respectively, N denotes the number of phase shift patterns, N ≥ 3, i ∈ {1, 2, ..., N}, a(x p ,y p ) is the average coding strength, 0<a(x p ,y p )≤255; b(x p ,y p ) is the coding modulation intensity, 0<b(x p ,y p )≤255. p ,y p ) and b(x p ,y p ) is 127.5. ]

[0024] Step 2.1, a projector projects a non-uniform periodic phase-shifted fringe pattern into the scene to be tested, and an industrial camera synchronously photographs the pattern to obtain a gray deformed fringe pattern, and further calculates the nonlinear wrapping phase;

[0025] The gray deformed fringe diagram is represented as:

[0026]

Formula derivation process: α(x c ,y c ) represents the reflectivity of the object. c ,y c ) represents the intensity of ambient light on the surface of the object, β2(x c ,y c ) represents the intensity of ambient light directly entering the camera from the environment. Indicates the gray distorted striped pattern captured by the camera.

[0027] Among them, x c and c Respectively represent the horizontal and vertical coordinates of the camera plane; A(x c ,y c ) is the average pixel intensity of the entire pattern set; A(x c ,y c ) can be represented by A. The image coordinates are only used to maintain the unity of the pixel variables in the formula and have no other meaning. c ,y c ) represents the intensity modulation of a given pixel; φ DB (x c ,y c ) represents the nonlinear wrapping phase;

[0028] The nonlinear wrapping phase calculation formula is:

[0029] in, Represents a gray distorted striped pattern captured by a camera.

[0030] Step 3 includes:

[0031] Step 3.1, calculate the nonlinear wrapping phase along y c Nonlinear wrapped phase gradient and nonlinear absolute phase gradient in increasing direction;

[0032] Step 3.2, construct an approximate phase using the nonlinear phase gradient, and use the approximate phase to unwrap the nonlinear local phase to obtain the nonlinear local phase;

[0033] The nonlinear wrapping phase gradient calculation formula in step 3.1 is:

[0034]

[0035] Specifically,

[0036] in, represents the gradient operator, Indicates along y c Nonlinear wrapped phase gradient in increasing direction;

[0037] The nonlinear absolute phase gradient calculation formula is:

[0038]

[0039] in, Indicates along y c Nonlinear absolute phase gradient in increasing direction;

[0040] The approximate phase calculation formula in step 3.2 is:

[0041]

[0042] in

[0043]

[0044] represents the approximate phase, ρ represents the approximate phase gradient, L>(n-1)×max{f1,f2,...f m}.

[0045] The nonlinear local phase calculation formula in step 3.2 is:

[0046]

[0047] in, represents the nonlinear local phase.

[0048] The linear absolute phase in step 4 is obtained as follows:

[0049] Step 4.1, obtaining a reliable phase order by minimizing the nonlinear error in the frequency domain of the candidate phase;

[0050] Step 4.2, the phase order is combined with the nonlinear wrapped phase to perform phase unwrapping to obtain the nonlinear absolute phase, and then the linear absolute phase is calculated through piecewise linear transformation;

[0051] The phase order calculation formula in step 4.1 is:

[0052]

[0053] in

[0054] F[κ]=∑<Fitting{Fd[κ]}-Fd[κ]>

[0055] Fd[κ]=log{abs[FT(κ)]}

[0056] k(x c ,y c ) represents the phase order, FT(.) represents the Fourier transform operator, abs[.] represents the absolute value operator, log{.} represents the logarithmic operator, and Fd[.] represents the transformation from the time domain to the frequency domain. F[.] represents the nonlinear error value, and: F[κ] = ∑<Fitting{Fd[κ]}-Fd[κ]> , Fitting{.} represents a two-dimensional nonlinear fitting operation, argmin{.} represents a minimum value search operator; the nonlinear absolute phase calculation formula in step 4.2 is:

[0057] Φ DB (x c ,y c )=φ DB (x c ,y c )+2π×k(x c ,y c )

[0058] Among them, Φ DB (x c ,y c ) represents the nonlinear absolute phase.

[0059] The linear absolute phase calculation formula in step 4.2 is:

[0060] Φ(x c ,y c )=H(Φ DB (x c ,y c))

[0061] in

[0062]

[0063] Among them, Φ(x c ,y c ) represents a linear absolute phase with a constant gradient, floor[.] represents a floor operator, and mod[.] represents a modulo operator. db_i ∈F DB , f p Indicates the fringe frequency used during the projector calibration phase.

[0064] Beneficial effects:

[0065] 1. The present invention proposes a phase unwrapping method for non-equal period phase shift, which uses multiple single-frequency fringes with different frequencies to perform combined encoding according to the De Bruijn order. Our method uses only 1 / 3 of the number of projections to achieve measurement accuracy similar to that of the advanced multi-frequency heterodyne method, greatly improving the coding efficiency;

[0066] 2. The present invention selects m kinds of sinusoidal fringes with similar wavelengths for combined coding, which ensures the uniformity of coding accuracy. Pure sinusoidal coding is used, which is robust to complex structures and different reflection areas of objects.

[0067] 3. The present invention does not require the introduction of additional equipment and does not require the projection of additional stripes, thereby ensuring low cost and high measurement efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 It is a flow chart of the phase unwrapping method of non-uniform periodic phase shift in Embodiment 1 of the present invention;

[0069] Figure 2 The three-step De Bruijn phase-shift fringe coding pattern in Embodiment 1 of the present invention;

[0070] Figure 3 These are three gray deformed fringe images obtained in Example 1 of the present invention;

[0071] Figure 4 is the nonlinear wrapped phase diagram in Embodiment 1 of the present invention;

[0072] Figure 5 is the nonlinear wrapped phase gradient image in Example 1 of the present invention;

[0073] Figure 6 is the nonlinear absolute phase gradient map in Example 1 of the present invention;

[0074] Figure 7It is a cross-sectional diagram of the approximate phase and the nonlinear local phase in Example 1 of the present invention.

[0075] Figure 8 This is the linear absolute phase diagram in Example 1 of the present invention.

[0076] Fig. 9 This is the object model in Example 1 of the present invention. DETAILED DESCRIPTION

[0077] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments. However, this should not be understood as the scope of the above subject matter of the present invention being limited to the following embodiments, and all technologies realized based on the content of the present invention belong to the scope of the present invention.

[0078] In order to better illustrate the method of the present invention, this embodiment takes a 5-element 3rd-order De Bruijn sequence as an example, combined with a three-step phase shift as an encoding scheme, and illustrates the measurement of a complex scene with a piece of aircraft skin part with rivets, a step block, a plastic bottle, two complex plaster figures and a pair of plaster dental models.

[0079] References Figure 1 The flowchart shown in the figure shows a phase unwrapping method for non-equal period phase shift, comprising the following steps:

[0080] Step 1: Construct a wavelength-De Bruijn element mapping table to generate a non-uniform periodic phase-shifted fringe group.

[0081] In order to ensure the uniformity of the calculated phase, five approximate wavelengths are selected for sinusoidal encoding, including: 12, 14, 16, 18 and 20. Then the wavelengths are assigned to five De Bruijn elements to construct a wavelength-De Bruijn element mapping table.

[0082] De Bruijn element-wavelength mapping table

[0083]

[0084]

[0085] Then, a set of 5-element 3rd-order De Bruijn wavelength sequences is constructed, as shown in formula (1);

[0086]

[0087] Among them, F DB represents a 5-ary 3rd order De Bruijn sequence. Next, use F DB The wavelength in generates a three-step De Bruijn phase-shifted fringe pattern, such as Figure 2 As shown, the calculation formula is:

[0088]

[0089] in

[0090]

[0091] x p and p Respectively represent the horizontal and vertical coordinates of the projector plane, I i DB (x p ,y p ) is the phase shift pattern, i∈{1,2,3}, and I i (x p ,y p ) and burned into the DLP projector's memory.

[0092] Step 2: Project the phase-shifted fringes frame by frame, obtain the modulation pattern and decode it to get the nonlinear wrapped phase.

[0093] The projector projects a three-step De Bruijn phase-shifted fringe pattern onto the scene to be tested, and the industrial camera simultaneously captures the gray deformed fringe pattern. like Figure 3 As shown. Among them, x c and c They represent the horizontal and vertical coordinates of the camera plane respectively, and the nonlinear wrapping phase is further calculated using formula (4);

[0094]

[0095] Among them, φ DB (x c ,y c ) represents the nonlinear wrapping phase, such as Figure 4 shown.

[0096] Step 3: Construct an approximate phase and locally expand the nonlinear wrapped phase to obtain the nonlinear local phase.

[0097] First, the nonlinear wrapping phase φ is calculated using formula (5) DB (x c ,y c ) along y c Gradient in increasing direction.

[0098]

[0099] in, represents the gradient operator, Indicates along y c The nonlinear wrapped phase gradient in increasing direction, such as Figure 5 As shown in Figure 6. Since the nonlinear wrapping phase has a phase truncation of 2π amplitude, the phase gradient will have periodic jumps. Therefore, we use formula (6) to Do gradient correction

[0100]

[0101] in, Indicates along y c The nonlinear absolute phase gradient in the increasing direction, such as Figure 6 Next, The local mean of is used as the approximate gradient to generate the approximate phase

[0102]

[0103] in

[0104]

[0105] ρ represents the approximate phase gradient, Represents the approximate phase, such as Figure 7 As shown. L>(n-1)×max{f1,f2,...f m}. Next, use the approximate phase To assist the nonlinear wrapping phase φ DB Do local phase unwrapping.

[0106]

[0107] in, represents the unwrapped nonlinear local phase, and The cross section of Figure 7 shown.

[0108] Step 4: Minimize the nonlinear error on the mapped phase to calculate the phase order, combine the nonlinear wrapped phase for phase unwrapping, and obtain the linear absolute phase through piecewise linear mapping.

[0109] Since only when With the correct k(x c ,y c ) can map the correct phase, otherwise there will be additional nonlinear errors on the phase surface. Therefore, each k(x c ,y c ) Linear phase mapping is performed to obtain a large number of candidate phases, and a reliable phase order is obtained by minimizing the nonlinear error in the frequency domain of the candidate phases. The calculation formula is:

[0110]

[0111] in

[0112] F[κ]=∑<Fitting{Fd[κ]}-Fd[κ]> (11)

[0113] Fd[κ]=log{abs[FT(κ)]} (12)

[0114] k(x c ,y c ) represents the phase order, FT(.) represents the Fourier transform operator, abs[.] represents the absolute value operator, log{.} represents the logarithmic operator, and Fd[.] represents the transformation from the time domain to the frequency domain. F[.] represents the nonlinear error value, Fitting{.} represents the two-dimensional nonlinear fitting operation, and argmin{.} represents the minimum value search operator. Next, the nonlinear wrapped phase φ DB (x c ,y c ) combined with the phase order k(x c ,y c )Use formula (13) for phase unwrapping.

[0115] Φ DB (x c ,y c )=φ DB (x c ,y c )+2π×k(x c ,y c ) (13)

[0116] Among them, Φ DB (x c ,y c ) represents the nonlinear absolute phase. Then, the nonlinear absolute phase is mapped into the linear absolute phase using the piecewise linear transformation of formula (14-16).

[0117] Φ(x c ,y c )=H(Φ DB (x c ,y c )) (14)

[0118] in

[0119]

[0120] Among them, Φ(x c ,y c ) represents a linear absolute phase with a constant gradient, such as Figure 8As shown. floor[.] represents the floor operator, and mod[.] represents the remainder operator. db_i ∈F DB , f p represents the fringe frequency used in the projector calibration stage. Finally, the linear absolute phase is combined with the structured light system parameters to reconstruct the object model (see reference: Li, Beiwen, et al. "High-accuracy, high-speed 3D structured light imaging techniques and potential applications to intelligent robotics." International journal of intelligent robotics and applications 1.1 (2017): 86-103.), such as Fig. 9 shown.

Claims

1. A three-dimensional measurement method of non-uniform periodic phase shift, characterized in that: The following steps are involved: Step 1: Construct a wavelength-De Bruijn element mapping table to generate a non-uniform periodic phase-shifted fringe group; Step 2: Project the phase-shifted fringes onto the object under test frame by frame, obtain the modulation pattern and decode it to obtain the nonlinear wrapped phase; Step 3: Construct an approximate phase and perform a local phase unwrapping operation on the nonlinear wrapped phase to obtain a nonlinear local phase; Step 4: Calculate the phase order by minimizing the nonlinear error on the nonlinear local phase, then combine the phase order with the nonlinear wrapped phase for phase unfolding, and then obtain the linear absolute phase of the object under test through piecewise linear mapping to complete the three-dimensional measurement.

2. The three-dimensional measurement method of non-uniform periodic phase shift according to claim 1, characterized in that: Step 1 includes: Step 1.1, select m wavelengths and De Bruijn elements, and construct a wavelength-De Bruijn element mapping table; The wavelength-De Bruijn element mapping table is: Step 1.2, generating an m-element n-order De Bruijn wavelength sequence and generating a non-uniform periodic phase shift code; The De Bruijn wavelength sequence is: The non-equal period phase shift encoding formula in step 1.2 is: in is the phase shift pattern, x p and p denote the horizontal and vertical coordinates of the projector plane, respectively, N denotes the number of phase shift patterns, N ≥ 3, i ∈ {1, 2, ..., N}, a(x p ,y p ) is the average coding strength, 0<a(x p ,y p )≤255; b(x p ,y p ) is the coding modulation intensity, 0<b(x p ,y p )≤255. p ,y p ) and b(x p ,y p ) is 127.

5. ] Step 2 includes: Step 2.1, a projector projects a non-uniform periodic phase-shifted fringe pattern into the scene to be tested, and an industrial camera synchronously photographs the pattern to obtain a gray deformed fringe pattern, and further calculates the nonlinear wrapping phase; The gray deformed fringe diagram is represented as: Among them, x c and c Respectively represent the horizontal and vertical coordinates of the camera plane; A(x c ,y c ) is the average pixel intensity B(x) of the entire pattern set c ,y c ) represents the intensity modulation of a given pixel; φ DB (x c ,y c ) represents the nonlinear wrapping phase; The nonlinear wrapping phase calculation formula is: in, Represents a gray distorted striped pattern captured by a camera. Step 3 includes: Step 3.1, calculate the nonlinear wrapping phase along y c Nonlinear wrapped phase gradient and nonlinear absolute phase gradient in increasing direction; Step 3.2, construct an approximate phase using the nonlinear phase gradient, and use the approximate phase to unwrap the nonlinear local phase to obtain the nonlinear local phase; The nonlinear wrapping phase gradient calculation formula in step 3.1 is: in, represents the gradient operator, Indicates along y c Nonlinear wrapped phase gradient in increasing direction; The nonlinear absolute phase gradient calculation formula is: in, Indicates along y c Nonlinear absolute phase gradient in increasing direction; The approximate phase calculation formula in step 3.2 is: in represents the approximate phase, ρ represents the approximate phase gradient, L>(n-1)×max{f1,f2,...f m }; The nonlinear local phase calculation formula in step 3.2 is: in, represents the nonlinear local phase. The linear absolute phase in step 4 is obtained as follows: Step 4.1, obtaining a reliable phase order by minimizing the nonlinear error in the frequency domain of the candidate phase; Step 4.2, the phase order is combined with the nonlinear wrapped phase to perform phase unwrapping to obtain the nonlinear absolute phase, and then the linear absolute phase is calculated through piecewise linear transformation; The phase order calculation formula in step 4.1 is: in F[κ]=∑<Fitting{Fd[κ]}-Fd[κ]> Fd[κ]=log{abs[FT(κ)]} k(x c ,y c ) represents the phase order, FT(.) represents the Fourier transform operator, abs[.] represents the absolute value operator, log{.} represents the logarithmic operator, Fd[.] represents the transformation from time domain to frequency domain. F[.] represents the nonlinear error value, Fitting{.} represents the two-dimensional nonlinear fitting operation, and argmin{.} represents the minimum value retrieval operator; The nonlinear absolute phase calculation formula in step 4.2 is: Φ DB (x c ,and c )=φ DB (x c ,and c )+2π×k(x c ,and c ) Among them, Φ DB (x c ,y c ) represents the nonlinear absolute phase. The linear absolute phase calculation formula in step 4.2 is: Φ(x c ,y c )=H(Φ DB (x c ,y c )) in Among them, Φ(x c ,y c ) represents a linear absolute phase with a constant gradient, floor[.] represents a floor operator, and mod[.] represents a modulo operator; f db_i ∈F DB , f p Indicates the fringe frequency used during the projector calibration phase.

Images